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Introduction 


Network  flow  problems  are  central  problems  in  operations  research,  computer  science,  and  engineer¬ 
ing  and  they  arise  in  many  real  world  applications.  Starting  with  early  work  in  linear  programming 
and  spurred  by  the  classic  book  of  Ford  and  Fulkerson^f^fl]',  the  study  of  such  problems  has  led  to 
continuing  improvements  in  the  efficiency  of  network  flow  algorithms.  In  spite  of  the  long  history 
of  this  study,  many  substantial  results  have  been  obtained  within  the  last  several  years.  In  this 
survey  we  examine  some  of  these  recent  developments  and  the  ideas  behind  them. 

We  discuss  the  classical  network  flow  problems,  the  maximum  flow  problem  and  the  rninimum- 
cost  circulation  problem ,  and  a  less  standard  problem,  the  generalized  flow  problem ,  sometimes 
called  the  problem  of  flows  with  losses  and  gains.  The  survey  contains  six  chapters  in  addition  to 
this  introduction.  Chapter  1  develops  the  terminology  needed  to  discuss  network  flow  problems. 
Chapter  2  discusses  the  maximum  flow  problem.  Chapters  3,  4,  and  5  discuss  different  aspects 
of  the  minimum-cost  circulation  problem,  and  Chapter  6  discusses  the  generalized  flow  problem. 
In  the  remainder  of  this  introduction,  we  mention  some  of  the  history  of  network  flow  research, 
comment  on  some  of  the  results  to  be  presented  in  detail  in  later  sections,  and  mention  some  results 
not  covered. in  this  survey. 

We  axe  interested  in  algorithms  whose  running  time  is  small  as  a  function  of  the  size  of  the 
network  and  the  numbers  involved  (e.g.  capacities,  costs,  or  gains).  As  a  measure  of  the  network 
size,  we  use  n  to  denote  the  number  of  vertices  and  m  to  denote  the  number  of  arcs.  As  measures 
of  the  number  sizes,  we  use  U  to  denote  the  maximum  arc  capacity,  C  to  denote  the  maximum  arc 
cost,  and  B  (in  the  case  of  the  generalized  flow  problem)  to  denote  the  maximum  numerator  or 
denominator  of  the  arc  capacities  and  gains.  In  bounds  using  U,  C ,  or  B,  we  make  the  assumption 
that  the  capacities  and  costs  are  integers,  and  that  the  gains  in  the  generalized  flow  problem  are 
rational  numbers. 

We  are  most  interested  in  polynomial-time  algorithms.  We  make  the  following  distinctions. 
An  algorithm  is  pseudopolynomial  if  its  running  time  has  a  bound  that  is  polynomial  in  n,m . 
and  the  appropriate  subset  of  U,  C,  and  B.  An  algorithm  is  polynomial  if  its  running  time  has 
a  bound  that  is  polynomial  in  n,m,  and  the  appropriate  subset  of  log  U,  log  C,  and  log/?1.  An 

J  Ail  logarithms  in  this  paper  without  an  explicit  base  are  base  two. 
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algorithm  is  strongly  polynomial  if  its  running  tinm  has  a  bound  that  is  polynomial  in  n  and  m2. 
When  comparing  polynomial  algorithms  to  strongly  polynomial  ones  we  shall  use  the  similarity 
assumption  that  log  U,  log  C  and  logf?  are  0(logn)  [32].  We  shall  also  be  interested  in  strongly 
polynomial  algorithms,  however. 

The  network  flow  problems  discussed  in  this  survey  are  special  cases  of  linear  programming, 
and  thus  they  can  be  solved  by  general-purpose  linear  programming  algorithms.  However,  the 
combinatorial  structure  of  these  problems  makes  it  possible  to  obtain  more  efficient  algorithms. 

We  shall  not  discuss  in  detail  algorithms  that  are  based  on  general  linear  programming  methods. 
We  should  mention,  however,  that  the  first  algorithm  designed  for  network  flow  problems  was  the 
network  simplex  method  of  Dantzig  [19].  It  is  a  variant  of  the  Unear  programming  simplex  method 
designed  to  take  advantage  of  the  combinatorial  structure  of  network  flow  problems.  Variants  of 
the  simplex  method  that  avoid  cycling  give  an  exponential  bound  on  the  complexity  of  all  the 
network  flow  problems.  (Cunningham  [18]  gives  an  elegant  anti-cycling  strategy  for  the  network 
simplex  method  based  on  graph-theoretic  properties  of  the  minimum-cost  circulation  problem). 
Recently,  Goldfarb  and  Hao  [51]  have  designed  a  variant  of  the  primal  network  simplex  method 
for  the  maximum  flow  problem  that  runs  in  strongly  polynomial  time  (see  Section  2.1).  Orlin  [80] 
designed  a  variant  of  the  dual  network  simplex  method  for  the  minimum-cost  circulation  problem 
that  runs  in  strongly  polynomial  time  (see  Chapter  5).  For  along  time,  the  network  simplex  method 
has  been  the  method  of  choice  in  practice,  in  particular  for  the  minimum- »_ost  circulation  problem 
(see  e.g.  [53]);  for  large  instances  of  hard  problems,  the  new  scaling  algorithms  are  probably  better, 
however. 

The  first  pseudopolynomial  algorithm  for  the  maximum  flow  problem  is  the  augmenting  path 
algorithm  of  Ford  and  Fulkerson  [25,  26].  Dinic  [20]  and  Edmonds  and  Karp  [21]  independently 
obtained  polynomial  versions  of  the  augmenting  path  algorithm.  Since  then,  several  more-efficient 
algorithms  have  been  developed.  Chapter  2  presents  the  push/relabel  method ,  recently  proposed  by 
Goldberg  [39]  and  Goldberg  and  Tarjan  [44],  along  with  some  of  its  more  efficient  variants. 

The  first  pseudopolynomial  algorithm  for  the  minimum-cost  circulation  problem  is  the  out- 
of-kilter  method,  which  was  developed  independently  by  Yakovleva  [104],  Minty  [76],  and  Fulk¬ 
erson  [31].  The  first  polynomial  algorithm  for  the  minimum-cost  circulation  problem  is  due  to 
Edmonds  and  Karp  [21].  To  develop  this  algorithm  Edmonds  and  Karp  introduced  the  technique 
of  scaling ,  which  has  proved  to  be  a  useful  tool  in  the  design  and  analysis  of  algorithms  for  a 
variety  of  combinatorial  optimization  problems.  Chapter  3  and  Section  5.2  are  devoted  to  scaling 
algorithms  for  the  minimum-cost  circulation  problem. 

The  maximum  flow  algorithms  of  Dinic  [20]  and  Edmonds  and  Karp  [21]  are  strongly  poly¬ 
nomial,  but  the  minimum-cost  circulation  algorithm  of  Edmonds  and  Karp  [21]  is  not.  The  first 
strongly  polynomial  algorithm  for  the  minimum-cost  circulation  problem  was  designed  by  Tar- 

2For  a  more  forma]  definition  of  polynomial  and  strongly  polynomial  algorithms,  see  [54]. 
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dos  [95].  Chapter  4  and  Section  5.3  are  devoted  to  recent  strongly  polynomial  algorithms  for  the 
minimum-cost  circulation  problem. 

The  first  augmenting  path  algorithms  for  the  generalized  flow  problem  were  developed  inde¬ 
pendently  by  Jewell  [61,  62]  and  Onaga  [78].  Many  pseudopolynomial  minimum-cost  circulation 
algorithms  have  been  adapted  for  the  generalized  flow  problem  (see  [101]  for  a  survey).  The  first 
polynomial-time  algorithm  for  the  generalized  flow  problem  was  the  ellipsoid  method  [69].  Kapoor 
and  Vaidya  [64]  have  shown  how  to  speed  up  Karmarkar  [65]  —  or  Renegar  [88]  —  type  interior- 
point  algorithms  on  network  flow  problems  by  taking  advantage  of  the  special  structure  of  the 
matrices  used  in  the  linear  programming  formulations  of  these  problems.  Vaidya’s  algorithm  [102] 
is  the  fastest  currently  known  algorithm  for  the  generalized  flow  problem.  The  first  polynomial 
algorithms  for  the  generalized  flow  problem  that  are  not  based  on  general-purpose  linear  program¬ 
ming  methods  are  due  to  Goldberg,  Plotkin,  and  Tardos  [42].  These  algorithms  are  discussed  in 
Chapter  6.  The  existence  of  a  strongly  polynomial  algorithm  for  the  generalized  flow  problem  is 
an  interesting  open  question. 

Important  special  cases  of  network  flow  problems  that  will  not  be  covered  in  this  survey  are  the 
bipartite  matching  problem  and  its  weighted  version,  the  assignment  problem.  These  problems  can 
be  stated  as  maximum  flow  and  minimum-cost  circulation  problems,  respectively,  on  networks  with 
unit  capacities  and  a  special  structure  (see  e.g.  [23,  96]).  Some  of  the  efficient  algorithms  for  the 
more  general  problems  have  evolved  from  efficient  algorithms  developed  earlier  for  these  simpler 
problems. 

Konig’s  [71]  proof  of  a  good  characterization  of  the  maximum  size  of  a  matching  in  a  bipartite 
graph  gives  an  0(nm)-time  algorithm  for  finding  a  maximum  matching.  The  Ford- Fulkerson 
maximum  flow  algorithm  can  be  viewed  as  an  extension  of  this  algorithm.  Hopcroft  and  Karp  [57] 
gave  an  0(\/nm )  algorithm  for  the  bipartite  matching  problem.  Even  and  Tarjan  observed  [24] 
that  Dinic’s  maximum  flow  algorithm,  when  applied  to  the  bipartite  matching  problem,  behaves 
similarly  to  the  Hopcroft-Karp  algorithm  and  runs  in  0(y/nm)  time  as  well.  A  variation  of  the 
Goldberg-Tarjan  maximum  flow  algorithm  (wdiich  can  be  viewed  as  a  generalization  of  Dinic's 
algorithm)  can  be  easily  shown  to  lead  to  the  same  bound  [5,  83].  In  spite  of  recent  progress  on 
related  problems,  the  O(yfnm)  bound  has  not  been  improved. 

The  first  algorithm  for  the  assignment  problem  is  the  Hungarian  method  of  Kuhn  [72].  The  out- 
of- kilter  algorithm  is  an  extension  of  this  algorithm  to  the  minimum-cost  circulation  problem.  The 
Hungarian  method  solves  the  assignment  problem  in  0(n )  shortest  path  computations.  Edmonds 
and  Karp  [21]  and  Tomizawa  [100]  have  observed  that  the  dual  variables  can  be  maintained  so  that 
these  shortest  path  computations  are  on  graphs  with  non-negative  arc  costs.  Combined  with  the 
shortest  path  algorithm  of  [29],  this  observation  gives  an  0(n(m  4-  n  log  n))  bound  for  the  problem. 
Gabow  [32]  used  scaling  to  obtain  an  C^n^m  log  C)  algorithm  for  the  problem.  Extending  ideas  of 
the  Hopcroft-Karp  bipartite  matching  algorithm  and  those  of  the  Goldberg-Tarjan  minimum-cost 
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Problem 

Date 

Discoverer 

Running  Time 

References 

Bipartite 

Matching 

1973 

Hopcroft  and  Karp 

Ofv/iim) 

[57] 

Assignment 

1955 

Kuhn 

0(n(m  +  n  log  n)) 

— 

1987 

Gabow  and  Tarjan 

0(y/n'l i  log(.iC)) 

mm 

Maximum 

1986 

Goldberg  and  Tarjan 

0(nm  lo g(n2/m)) 

[44] 

Flow 

1988 

Ahuja,  Orlin,  and 
Tarjan 

0(nm  log(%y/logU  +  2)) 

[4] 

Minimum-Cost 

1972 

Edmonds  and  Karp 

0(m  log  U(m  +  n  log  n)) 

[21] 

Circulation 

1987 

Goldberg  and  Tarjan 

0(nm  log (n2/m)  log(nC)) 

[46] 

1987 

Orlin 

0(m  log  n(m  +  n  log  n)) 

[82] 

1988 

Ahuja,  Goldberg, 
Orlin.  and  Tarjan 

0(nm  log  log  U  log(nC)) 

[1] 

Generalized 

Flow 

1989 

Vaidya 

0(n2ml  5  log(nB)) 

[102] 

Table  1:  Fastest  currently  known  algorithms  for  network  flow  problems. 


circulation  algorithm  (discussed  in  Section  3),  Gabow  and  Tarjan  obtained  an  0(^/nm\og{nC)) 
algorithm  for  the  assignment  problem. 

A  more  recent  pseudopolynomial  algorithm  for  the  assignment  problem  is  the  auction  algorithm 
of  Bertsekas  [9]  (first  published  in  [10]).  This  algorithm  contains  some  of  the  elements  of  the 
push/relabel  algorithms  discussed  in  Sections  2  and  3. 

Some  versions  of  the  network  simplex  method  have  been  shown  to  solve  the  assignment  problem 
in  polynomial  time.  In  particular,  Orlin  [81]  shows  that  a  natural  version  of  the  primal  simplex 
method  runs  in  polynomial  time,  and  Balinski  [6]  gives  a  signature  method  that  is  a  dual  simplex 
algorithm  for  the  assignment  problem  and  runs  in  strongly  polynomial  time. 

For  discussion  of  parallel  algorithms  for  the  bipartite  matching  and  assignment  problems,  see 

[10,  33,  43,  67,  77], 

In  this  survey  we  would  like  to  highlight  the  main  ideas  involved  in  designing  highly  efficient 
algorithms  for  network  flow  problems,  rather  than  to  discuss  in  detail  the  fastest  algorithms.  How¬ 
ever,  for  easy  reference,  we  summarize  the  running  times  of  the  fastest  currently  known  algorithms 
in  Table  1.  For  each  problem  we  list  all  of  the  bounds  that  are  best  for  some  values  of  the  param¬ 
eters,  but  we  only  list  the  first  algorithm  achieving  the  same  bound.  Some  of  the  bounds  stated  in 
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the  table  depend  on  the  0(m  +  nlogn)  implementation  of  Dijkstra’s  shortest  path  algorithm  [29]. 


Chapter  1 


Preliminaries 


In  this  chapter  we  define  the  problems  addressed  in  this  survey  and  review  fundamental  facts  about 
these  problems.  These  problems  are  the  maximum  flow  problem,  the  minimum-cost  circulation 
problem,  the  transshipment  problem,  and  the  generalized  flow  problem 


1.1  Flows  and  Residual  Graphs 

A  network  is  a  directed  graph  G  =  (V,E)  with  a  non-negative  capacity  function  u  :  E  — *  R.x.,* 
We  assume  that  G  has  no  multiple  arcs,  i.e.,  E  C  V  x  V .  If  there  is  an  arc  from  a  vertex  r  to  a 
vertex  w,  this  arc  is  unique  by  the  assumption,  and  we  will  denote  it  by  (v,  w).  This  assumption  is 
for  notational  convenience  only.  We  also  assume,  without  loss  of  generality,  that  the  input  graph 
G  is  symmetric:  (v,w)  £  E  <=>  (w,v)  £  E.  A  flow  network  is  a  network  with  two  distinguished 
vertices,  the  source  s  and  the  sink  t. 

A  pseudoflow  is  a  function  f  :  E  R  that  satisfies  the  following  constraints: 

f(v,w)  <  u(v,v>)  V(r, w)  £  E  (capacity  constraint),  (1.1) 

/( r,tr)  =  -  /(??, r)  V(r,u’)  £  E  (flm>-  antisymmetry  constraint).  (1.2) 

Remark:  To  gain  intuition,  it  is  often  useful  to  think  only  about  the  non-negative  components  of  a 
pseudoflow  (or  of  a  generalized  pseudoflow,  defined  below).  The  antisymmetry  constraints  reflect 
the  fact  that  a  flow  of  value  x  from  r  to  rr  can  be  thought  of  as  a  flow  of  value  (— i)  from  w  to  r. 
The  negative  flow  values  are  introduced  only  for  notational  convenience.  Note,  for  example,  that 
one  does  not  have  to  distinguish  between  lower  and  upper  capacity  bounds:  the  capacity  of  the  arc 
(v,  w)  represents  a  lower  bound  on  the  flow  value  on  the  opposite  arc. 

________  - 
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Given  a  pseud^r  »w  /,  we  define  the  excess  function  ej:V  — *  R  by  ej(v)  —  /(u,r),  the 

net  flow  into  v.  We  will  say  that  a  vertex  v  has  excess  if  e/(v)  is  positive,  and  has  deficit  if  it  is 
negative,  For  a  vertex  v,  we  define  the  conservation  constraint  by 

e/(u)  =  0  (flow  conservation  constraint).  (1.3) 

Given  a  pseudoflow  /,  the  residual  capacity  function  Uf  :  E  — *  R  is  defined  by  vj(v.w)  = 
u(u,  tr)  -  f(v,w).  The  residual  graph  with  respect  to  a  pseudoflow  /  is  given  by  Gj  -  (V.Ej). 
where  E j  =  {(r,  w)  £  E\uj(v,w)  >  0}. 


1.2  The  Maximum  Flow  Problem 

To  introduce  the  maximum  flow  problem,  we  need  the  following  definitions  in  addition  to  the 
definitions  of  the  previous  section.  Consider  a  flow  network  (G,u,s,t).  A  prcflou'  is  a  pseudoflow 
/  such  that  the  excess  function  is  non-negative  for  all  vertices  other  than  s  and  t.  A  flow  f  on  G 
is  a  pseudoflow  satisfying  the  conservation  constraints  for  all  vertices  except  s  and  /.  The  value  \f\ 
of  a  flow  /  is  the  net  flow  into  the  sink  e j(t).  A  maximum  flow  is  a  flow  of  maximum  value  (also 
called  an  optimal  flow).  The  maximum  flow  problem  is  that  of  finding  a  maximum  flow  in  a  given 
flow  network. 

Given  a  flow  /,  we  define  an  augmenting  path  to  be  a  source-to-sink  path  in  the  residual  graph. 
The  following  theorem,  due  to  Ford  and  Fulkerson,  gives  an  optimality  criterion  for  maximum  flows. 

Theorem  1.2.1  [2G]  A  flow  is  optimal  if  and  only  if  its  residual  graph  contains  no  augmenting  path 


1.3  The  Minimum-Cost  Circulation  Problem 

A  circulation  is  a  pseudoflow  with  zero  excess  at  every  vertex.  A  cost  function  is  a  real-valued 
function  on  arcs  c  :  E  —• ■  R.  Without  loss  of  generality,  we  assume  that  costs  are  antisymmetric: 

c(r,  u1)  =  -c(u\  u)  V(r, tr)  £  E  (cost  antisymmetry  constraint ).  (1.4) 

The  cost  of  a  circulation  /  is  given  by 

c(f)=  f(-,w)c{v,w). 


The  minimum-cost  circulation  problem  is  that  of  finding  a  minimum-cost  (optimal)  circulation  in 
an  input  network  (G,u,c). 


1.3.  THE  MINIMUM-COST  CIRCULATION  PROBLEM 


11 


We  have  assumed  that  the  capacities  are  non-negative.  This  assumption  is  no  loss  of  generality. 
Given  a  circulation  problem  with  some  negative  capacities  one  can  find  a  feasible  circulation  /  using 
one  invocation  of  a  maximum  flow  algorithm.  (See  [86],  problem  11(e),  p.  215.)  The  non-negative 
residual  capacities  u(v,w)  -  f(v,w)  then  define  an  equivalent  problem. 

Next  we  state  two  criteria  for  the  optimality  of  a  circulation.  Define  the  cost  of  a  cycle  T  to 
be  the  sum  of  the  costs  of  the  arcs  along  the  cycle  and  denote  it  by  c(T).  The  following  theorem 
states  the  first  optimality  criterion. 

Theorem  1.3.1  [15]  A  circulation  is  optimal  if  and  only  if  its  residual  graph  contains  no  negative-cost 
cycle. 

To  state  the  second  criterion,  we  need  the  notions  of  a  price  function  and  a  reduced  cost  function. 
A  price  function  is  a  vertex  labeling  p  :  V  — *  R.  The  reduced  cost  function  with  respect  to  a  price 
function  p  is  defined  by  cp(v,  w)  —  c(v,  w)+p(v)  —  p(w).  These  notions,  which  originate  in  the  theory 
of  linear  programming,  are  crucial  for  many  minimum-cost  flow  algorithms.  As  linear  programming 
dual  variables,  vertex  prices  have  a  natural  economic  interpretation  as  the  current  market  prices 
of  the  commodity.  We  can  interpret  the  reduced  cost  cp(i\u')  as  the  cost  of  buying  a  unit  of 
commodity  at  v,  transporting  it  to  it,  and  then  selling  it.  Due  to  the  conservation  constraints,  the 
reduced  costs  define  an  equivalent  problem. 

Theorem  1.3.2  [26]  A  circulation  /  is  optimal  for  the  minimum-cost  circulation  problem  ( G.u.c )  if 
and  only  if  it  is  optimal  for  the  problem  ( G,U,cp )  for  every  price  function  p. 

The  second  optimality  criterion  is  as  follows. 

Theorem  1.3.3  [26]  A  circulation  /  is  optimal  if  and  only  if  there  is  a  price  function  p  such  that,  for 
each  arc  (r,  w), 

cp(v,w)  <  0  =>  f(i\x)  =  u{ v,w)  (complementary  slackness  constraint).  (1-5) 

A  minimum-cost  circulation  problem  may  have  no  optimal  solution.  The  following  theorem 
characterizes  when  this  is  the  case. 

Theorem  1.3.4  There  exists  a  minimum-cost  circulation  if  and  only  if  the  input  network  contains  no 
negative-cost  cycle  consisting  of  arcs  with  infinite  capacity. 

Note  that  the  maximum  flow  problem  is  a  special  case  ol  the  minimum-cost  circulation  problem 
(see  e.g.  [26]).  To  see  this,  consider  a  flow  network  and  add  a  pair  of  new  arcs  (s.t)  and  {t.s) 
with  u(s,t)  —  0  and  u(t,s>)  =  oo.  Define  the  costs  of  all  original  arcs  to  he  zero,  and  define 
c(s,t)  =  -c(t,s)  =  1.  A  minimum-cost  circulation  in  the  resulting  network  restricted  to  the 
original  network  is  a  maximum  flow  in  the  original  network. 


12 


CHAPTER  1.  PRELIMINARIES 


1.4  The  Transshipment  Problem 

Tn  this  section  we  define  the  (uncapacitated)  transshipment  problem.  Although  this  problem  is 
equivalent  to  the  minimum-cost  circulation  problem,  some  algorithms  are  more  natural  in  the 
context  of  the  transshipment  problem.  For  a  more  extensive  discussion  of  this  problem  and  the 
closely  related  transportation  problem,  see  [73]. 

In  the  transshipment  problem,  all  arc  capacities  are  either  zero  or  infinity.  In  addition,  the 
input  to  the  problem  contains  a  demand  function  d  :  V  — ►  R  such  that  Ylvev  d{v)  =  0.  For  the 
transshipment  problem,  the  notion  of  the  excess  at  a  vertex  is  defined  as  follows: 

eAv)  =  f(w'v)  ~  <1’)-  (1-6) 

wev 


A  i  seudoflow  /  is  feasible  if  the  conservation  constraints  (1.3)  hold  for  all  vertices.  The  transship¬ 
ment  problem  is  that  cr  finding  a  minimum-cost  feasible  pseudoflow  in  an  input  network.  In  the 
special  case  of  integer  demands,  we  shall  use  D  to  denote  the  maximum  absolute  value  of  a  demand. 

Theorems  analogous  to  Theorem  1.3.1,  1.3.3  and  1.3.2.  hold  for  the  transshipment  problem, 
and  the  analog  of  Theorem  1.3.4  holds  for  a  transshipment  problem  that  has  a  feasible  pseudoflow. 

We  make  two  simplifying  assumptions  about  the  transshipment  problem.  First,  we  assume  that, 
for  the  zero  flow,  all  residual  arcs  have  non-negative  cost.  This  assun.ption  can  be  validated  by  first 
checking  whether  the  condition  of  Theorem  1.3.4  is  satisfied  (using  a  shortest  path  computation). 
In  the  case  of  a  transshipment  problem,  the  residual  graph  of  the  zero  flow'  consists  of  the  arcs 
w'ith  infinite  capacity.  If  this  graph  has  no  negative-cost  cycles,  then  define  p(v)  to  be  the  cost  of  a 
minimum-cost  path  from  a  fixed  vertex  s  to  the  vertex  v.  The  costs  cp(v ,w)  define  an  equivalent 
problem  that  satisfies  the  assumption.  Second,  we  assume  that  the  residual  graph  of  the  zero  flow  is 
strongly  connected.  (There  is  a  path  between  each  pair  of  vertices.)  This  condition  can  be  imposed 
by  adding  two  new  vertices  x  and  y  and  adding  an  appropriately  high-cost  arc  from  x  to  y  with 
u{x,y)  =  oo  and  u(y,x)  =  0,  and  arcs  between  the  new  vertices  and  every  other  vertex  in  both 
directions  such  that  u(v,x)  =  u(y,v)  =  oo  and  u(x,v)  =  u(v,y )  =  0  for  every  vertex  v.  If  the 
original  problem  has  a  feasible  solution,  then  in  every  minimum-cost  solution  of  the  transformed 
problem  all  of  the  dummy  arcs  introduced  to  make  the  graph  strongly  connected  have  zero  flow. 

Next  we  show'  that  the  transshipment  problem  is  equivalent  to  the  minimum-cost  circulation 
problem  (see  e.g.  [26]).  Given  an  instance  of  the  transshipment  problem,  we  construct  an  equivalent 
instance  of  the  minimum-cost  circulation  problem  as  follows.  Add  two  new  vertices,  x  and  y,  and 
arcs  (y,x)  and  (x,y)  with  u(x,y)  =  0  and  u (y,x)  =  J2v.d(v)>od(v).  Define  the  cost  of  (y,x)  to 
be  small  enough  so  that  any  simple  cycle  containing  ( y,x )  has  negative  cost;  for  example,  define 
c(y,x )  =  -c(x,y)  =  -nC.  For  every  vertex  v  €  V  with  d(v)  >  0,  add  arcs  (x,r)  and  (v.  x)  and 
define  u(x,v)  =  d(v),  u(v,x)  =  0,  and  c(x,v)  =  c(u,x)  =  0.  For  every  vertex  v  e  V  with  d(v)  <  0, 
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add  arcs  ( v,y )  and  (jt,v)  and  define  u(v,y)  =  -d(v),  u(y,v)  =  0,  and  c(v,y)  -  c{y,v)  =  0. 
Consider  an  optimal  solution  /  to  the  minimum-cost  circulation  problem  in  the  resulting  network. 
The  transshipment  problem  is  feasible  if  and  only  if  /  saturates  all  new  arcs,  and  in  this  case  the 
restriction  of  /  to  the  original  arcs  gives  an  optimal  solution  to  the  transshipment  problem. 

Next  we  reduce  the  minimum-cost  circulation  problem  to  the  transshipment  problem.  The  re¬ 
duction  uses  the  technique  of  edge-splitting.  Consider  a  minimum-cost  circulation  problem  (G,  u.  e). 
First  we  make  sure  that  for  every  edge  {u,  te}  of  the  network,  either  u(v,w)  or  u(w,v)  is  infinite. 
For  every  edge  {u,u’}  that  does  not  satisfy  the  above  condition,  we  introduce  a  new  vertex 
and  replace  the  arcs  (v,w)  and  (io,t;)  by  the  arcs  (z(vw),w),  (w,z(VtUl)),  and  (x(„,u.),  r). 

Define  the  costs  and  capacities  of  the  new  arcs  as  follows: 


u(v,Z(VlW))  =  u(v,w ), 
u(z(V'W),tr)  =  oo, 

tt(«\Z(v,u;))  = 
u(*(t ,,w)*v)  - 


c(v>X(v,w))  =  c(l\w), 
c(X(v,wbU’)  =  °’ 
C(»,l(v.u.))  =  0, 
C(*(v,u,)it>)  =  c(  U\  t). 


This  defines  an  equivalent  minimum-cost  circulation  problem  in  which  the  capacity  of  every  edge 
is  unbounded  in  at  least  one  direction. 

To  complete  the  transformation,  we  need  to  force  the  additional  restriction  that  every  finite 
capacity  is  zero.  Initialize  the  demand  function  d  to  be  zero  on  all  vertices.  Let  {i\tr}  be  an  edge 
such  that  u(w,v)  =  oo  and  u(v,w)  is  finite  but  not  zero.  Replacing  u(v,  it)  by  zero,  and  d(v) 
and  d(w)  by  d( v)  -f  u(v,w)  and  d(w)  —  u( v,w),  respectively,  defines  an  equivalent  problem  with 
u(v,  w)  =  0. 

Remark:  This  transformation  increases  the  number  of  vertices:  from  a  network  with  n  vertices 
and  m  edges,  we  obtain  a  network  with  up  to  n  +  m  vertices  and  2m  edges.  However,  the  resulting 
network  has  a  special  structure.  It  is  important  for  the  analyses  in  Chapter  5  to  notice  that  the 
newly  introduced  vertices  can  be  eliminated  for  shortest  path  computations.  Consequently,  this 
blowup  in  the  number  of  vertices  will  not  affect  the  running  time  of  shortest  path  computations  in 
the  residual  graph. 

The  capacitated  transshipment  problem  generalizes  both  the  minimum-cost  circulation  problem 
and  the  (uncapacitated)  transshipment  problem  discussed  above.  It  is  the  same  as  the  uncapac¬ 
itated  transshipment  problem  without  the  assumption  that  all  arc  capacities  are  infinite  or  zero. 
The  reductions  above  can  be  extended  to  show  that  the  capacitated  transshipment  problem  is 
equivalent  to  the  minimum-cost  circulation  problem.  Whereas  the  previous  simpler  formulations 
are  better  suited  for  designing  algorithms,  the  more  general  form  can  be  useful  in  applications. 
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In  this  section  we  define  generalized  pseudoflows,  generalized  flows,  and  the  generalized  flow  prob¬ 
lem,  and  compare  these  concepts  with  their  counterparts  discussed  in  the  previous  sections.  For 
alternative  formulations  of  the  problem,  see  e.g.  [42,  73).  (This  problem  is  also  known  as  the 
problem  of  flows  with  losses  and  gains.) 

The  generalized  flow  problem  is  given  by  a  network  with  a  distinguished  vertex,  the  source  s, 
and  a  gain  function  7  :  E  — ►  R+2  on  the  arcs.  We  assume  (without  loss  of  generality)  that  the  gain 
function  is  antisymmetric: 

-y(v,  w)  =  V(v,w)  £  E  (gain  antisymmetry  constraints).  (1.7) 

In  the  case  of  ordinary  flows,  if  f(v,w)  units  of  flow  are  shipped  from  v  to  w ,  f(i\w)  units 
arrive  at  u\  In  the  case  of  generalized  flows,  if  g(v,w)  units  of  flow  are  shipped  from  v  to  u\ 
7 (v,w)g(v,w)  units  arrive  at  ir.  A  generalized  pseudoflow  is  a  function  g  \  E  —*  R  that  satisfies 
the  capacity  constraints  (1.1)  and  the  generalized  antisymmetry  constraints: 

g{v,w)  =  v)g(w.v)  V(v,w)  £  E  (generalized  antisymmetry  constraints).  (1.8) 

The  gain  of  a  path  (cycle)  is  the  product  of  the  gains  of  the  arcs  on  the  path  (cycle).  For  a  given 
generalized  pseudoflow  g,  the  residual  capacity  and  the  residual  graph  Gg  =  (V,Eg)  are  defined  in 
the  same  way  as  for  pseudoflows.  The  excess  eg(v)  of  a  generalized  pseudoflow  g  at  a  vertex  r  is 
defined  by 

e5(u)  =  -  5Z  9{v,w).  (1.9) 

We  will  say  that  a  vertex  v  has  excess  if  ep(u)  is  positive  and  deficit  if  it  is  negative. 

A  generalized  flow  is  a  generalized  pseudoflow  that  satisfies  the  conservation  constraints  (1.3) 
at  all  vertices  except  at  the  source.  The  value  of  a  generalized  pseudoflow  g  is  the  excess  eff(s).  The 
generalized  flow  problem  is  that  of  finding  a  generalized  flow  of  maximum  value  in  a  given  network. 
In  our  discussion  of  this  problem,  we  shall  assume  that  the  capacities  are  integers,  and  that  each 
gain  is  a  rational  number  given  as  the  ratio  of  two  integers.  Let  B  denote  the  largest  integer  used 
in  the  specification  of  the  capacities  and  gains. 

A  flow-generating  cycle  is  a  cycle  whose  gain  is  greater  than  1,  and  a  flow-absorbing  cycle  is  a 
cycle  whose  gain  is  less  than  1.  Observe  that  if  one  unit  of  flow  leaves  a  vertex  i>  and  travels  along 
a  flow-generating  cycle,  more  than  one  unit  of  flow  arrives  at  v.  Thus  we  can  augment  the  flow 

JR+  denotes  the  set  of  positive  reals. 
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along  this  cycle;  that  is,  we  can  increase  the  excess  at  any  vertex  of  the  cycle  while  preserving  the 
excesses  at  the  other  vertices,  by  increasing  the  flow  along  the  arcs  in  the  cycle  (and  correspondingly 
decreasing  the  flow  on  the  opposite  arcs  to  preserve  the  generalized  antisymmetry  constraints). 

The  generalized  flow  problem  has  the  following  interpretation  in  financial  analysis.  The  com¬ 
modity  being  moved  is  money,  nodes  correspond  to  different  currencies  and  securities,  arcs  cor¬ 
respond  to  possible  transactions,  and  gain  factors  represent  the  prices  or  the  exchange  rates  (see 
Section  6.1).  From  the  investor’s  point  of  view,  a  residual  flow-generating  cycle  is  an  opportunity 
to  make  a  profit.  It  is  possible  to  take  advantage  of  this  opportunity,  however,  only  if  there  is  a 
way  to  transfer  the  profit  to  the  investor’s  bank  account  (the  source  vertex).  This  motivates  the 
following  definition.  A  generalized  augmenting  path  (GAP)  is  a  residual  flow-generating  cycle  and 
a  (possibly  trivial)  residual  path  from  a  vertex  on  the  cycle  to  the  source.  Given  a  generalized  flow 
and  a  GAP  in  the  residual  graph,  we  can  augment  the  flow  along  the  GAP,  increasing  the  value 
of  the  current  flow.  The  role  of  GAP’s  in  the  generalized  flow  problem  is  similar  to  the  role  of 
negative-cost  cycles  in  the  minimum-cost  circulation  problem  -  both  can  be  used  to  augment  the 
flow  and  improve  the  value  of  the  current  solution.  Onaga  [79]  proved  that  the  non-existence  of 
GAP’s  in  the  residual  graph  characterizes  the  optimality  of  a  generalized  flow. 

Theorem  1.5.1  [79]  A  generalized  flow  is  optimal  if  and  only  if  its  residual  graph  contains  no  GAP. 

Using  the  linear  programming  dual  of  the  problem,  it  is  possible  to  give  an  alternate  optimality 
criterion,  similar  to  the  criterion  in  Theorem  1.3.3  for  the  minimum-cost  circulation  problem.  A 
price  function  p  is  a  labeling  of  the  vertices  by  real  numbers,  such  that  p(s)  =  1.  As  in  the  case 
of  the  minimum-cost  circulation  problem,  vertex  prices  can  be  interpreted  as  market  prices  of  the 
commodity  at  vertices.  If  a  unit  of  flow  is  purchased  at  v  and  shipped  to  u\  then  7(1 \w)  units 
arrive  at  w.  Buying  the  unit  at  v  costs  p(v),  and  selling  7 [v,w)  units  at  w  returns  p(u')~){v,U'). 
Thus  the  reduced  cost  of  {v,w)  is  defined  as 

cp(v,w)  =  p(v )  -  p(w) -y(v,w). 

Linear  programming  duality  theory  provides  the  following  optimality  criterion. 

Theorem  1.5.2  [61]  A  generalized  flow  is  optimal  if  and  only  if  there  exists  a  price  function  p  such 
that  the  complementary  slackness  conditions  (1.5)  hold  for  each  arc  (v,w)  G  E. 

One  generalization  of  the  problem,  namely  the  generalized  flow  problem  with  costs ,  is  worth 
mentioning  here.  As  the  name  suggests,  in  this  problem  each  arc  (u,u>)  has  a  cost  c(v,w)  in 
addition  to  its  gain.  The  goal  is  to  find  a  generalized  flow  that  maximizes  eg(s)  -  c(g),  where  c{g)  = 
Yl(v,w):g(v,w)>oc(v'w)9(v'w)-  The  introduction  of  costs  enriches  the  combinatorial  structure  of  the 
problem  and  allows  the  modeling  of  more  complex  problems,  in  particular  economic  processes.  For 
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example,  a  positive  cost  flow-generating  cycle  with  a  path  leading  to  a  negative  cost  flow-absorbing 
cycle  may  be  used  as  a  producer-consumer  model.  The  generalized  flow  problem  with  costs  has  not 
been  studied  as  much  as  the  other  problems  discussed  in  this  survey,  and  the  only  polynomial-time 
algorithms  known  for  this  problem  are  based  on  general  linear  programming  methods  [64,  102]. 


1.6  The  Restricted  Problem 

Next  we  introduce  a  special  variant  of  the  generalized  flow  problem,  and  show  that  this  variant  is 
equivalent  to  the  original  problem. 

Consider  for  a  moment  the  following  variation  of  the  generalized  flow  problem:  given  a  flow 
network  w'ith  a  source  s  £  V  and  a  sink  t  £  V ,  find  a  generalized  pseudoflow  with  maxiir  :rr.  recess 
at  t  and  zero  excess  at  each  vertex  other  than  s  and  t.  Onaga  [78]  suggested  the  study  of  the  special 
case  of  this  problem  in  which  the  residual  graph  of  the  zero  flow  has  no  flow-generating  cycles.  We 
shall  consider  the  corresponding  special  case  of  the  generalized  flow  problem  in  which  the  residual 
graph  of  the  zero  flow  has  the  property  that  all  flow-generating  cycles  pass  through  the  source. 
(If  there  are  no  flow-generating  cycles,  the  zero  flow  is  the  optimal.)  We  shall  also  assume  that 
the  residual  graph  of  the  zero  flow  is  strongly  connected.  A  generalized  flow  network  in  which  the 
residual  graph  of  the  zero  flow  is  strongly  connected  and  which  has  no  flow-generating  cycles  not 
passing  through  the  source  is  called  a  restricted  network.  The  restricted  problem  is  the  generalized 
flow  problem  on  a  restricted  network.  The  restricted  problem  has  a  simpler  combinatorial  structure 
that  leads  to  simpler  algorithms.  Moreover,  it  turns  out  that  the  restricted  problem  is  equivalent 
to  the  generalized  flow  problem.  All  of  the  algorithms  that  we  review  solve  the  restricted  problem. 
In  the  rest  of  this  section  we  shall  review  basic  properties  of  the  restricted  problem,  and  we  outline 
the  reduction.  In  Chapter  6  we  will  use  the  term  “generalized  flow  problem”  to  mean  the  restricted 
problem. 

One  of  the  nice  facts  about  the  restricted  problem  is  that  the  optimality  condition  given  by 
Theorem  1.5.1  simplifies  in  this  case,  and  becomes  very  similar  to  Theorem  1.3.1.  This  character¬ 
ization,  which  is  also  due  to  Onaga  [78],  can  be  deduced  from  Theorem  1.5.1  with  the  use  of  the 
following  lemma.  The  lemma  can  be  proved  by  induction  on  the  number  of  axes  with  positive  flow. 

Lemma  1.6.1  Let  g  be  a  generalized  pseudoflow  in  a  restricted  network.  If  the  excess  at  every  vertex 
is  non-negative,  then  for  every  vertex  r  there  exists  a  path  from  v  to  s  in  the  residual  graph  Gg. 

Theorem  1.6.2  [79]  A  generalized  flow  g  in  a  restricted  network  is  optimal  if  and  only  if  the  residual 
graph  of  g  contains  no  flow-generating  cycles. 

Note  that  the  condition  in  the  theorem  can  be  formulated  equivalently  as  “if  and  only  if  the 
residual  graph  of  g  has  no  negative-cost  cycles,  where  the  cost  function  is  given  by  c  =  -  log 7.” 
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Theorem  1.6.3  [42]  The  generalized  flow  problem  can  he  reduced  to  the  restricted  problem  in  O(nm) 
time. 

Proof  sketch:  Given  an  instance  of  the  generalized  flow  problem,  reduce  it  to  an  instance  in  which 
all  the  flow-generating  cycles  pass  through  s,  as  follows.  Let  h  be  the  generalized  pseudoflow  that 
saturates  all  arcs  of  gain  greater  than  1  and  is  zero  on  all  arcs  of  gain  1.  Define  a  new  generalized 
flow  network  containing  each  arc  ( v ,  tn)  of  the  original  network,  with  new  capacity  u(v,  w)  -  h(v ,  w). 
For  each  vertex  v  £  V  such  that  e^(v)  >  0,  add  an  arc  (s, v)  of  capacity  u(s,v)  =  e/,(t’)/a,  and 
with  gain  7(5,  v)  =  a,  where  0  =  Bn.  Also  add  a  reverse  arc  (u,s)  with  u(r,s)  =  0.  For  each 
vertex  v  €  V  with  e/,(u)  <  0,  add  an  arc  (v,s)  of  capacity  u(u,s)  =  e^(t>)  and  with  gain  7 (v,s)  =  a  ; 
also  add  a  reverse  arc  (s,v)  with  u(v,  s)  =  0.  Let  g  be  an  optimal  generalized  flow  in  the  new 
network.  Consider  the  generalized  pseudoflow  g  +  h  restricted  to  the  arcs  of  the  original  network. 
Using  the  Decomposition  Theorem  1.7.3  below,  one  can  show  that  the  value  of  this  generalized 
pseudoflow  is  equal  to  the  value  of  the  optimal  generalized  flow,  and  an  optimal  generalized  flow 
can  be  constructed  from  this  generalized  pseudoflow  in  O(nm)  time.  Intuitively,  the  new  arcs  ensure 
that  vertices  having  excesses  with  respect  to  h  are  supplied  with  an  adequate  amount  of  “almost 
free”  flow,  and  vertices  having  deficits  with  respect  to  h  can  deliver  the  corresponding  amount  of 
flow  very  efficiently  to  the  source. 

Having  eliminated  flow-generating  cycles  not  containing  the  source,  we  can  impose  the  strong 
connectivity  condition  by  deleting  all  of  the  graph  except  the  strongly  connected  component  con¬ 
taining  the  source.  This  does  not  affect  the  value  of  an  optimum  solution  by  Theorem  1.5.1.  | 

Note  that  this  transformation  increases  the  size  of  the  numbers  in  the  problem  description. 
However,  the  polynomial-time  algorithms  described  later  depend  only  on  T  (<  Bn),  the  maximum 
product  of  gains  along  a  simple  path,  and  L  (<  B2m),  the  product  of  denominators  of  arc  gains; 
these  parameters  do  not  increase  significantly. 

1.7  Decomposition  Theorems 

A  useful  property  of  flows  and  generalized  flows  is  the  fact  that  they  can  be  decomposed  into  a 
small  number  of  “primitive  elements”.  These  elements  depend  on  the  problem  under  consideration. 
In  this  section  we  review  decomposition  theorems  for  circulations,  pseudoflows,  and  generalized 
pseudoflows. 

In  the  case  of  circulations,  the  primitive  elements  of  the  decomposition  are  flows  along  simple 
cycles. 

Theorem  1.7.1  [26]  For  every  circulation  /,  there  exists  a  collection  of  k  <  m  circulations  g\ , ...  ,<7*; 
such  that  for  every  i,  the  graph  induced  by  the  set  of  arcs  on  which  g,  is  positive  is  a  simple  cycle 
consisting  of  arcs  (r,  in)  with  f(v,w)  >  0,  and 
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Such  a  decomposition  can  be  found  in  O(nm)  time. 

For  pseudoflows  (or  flows),  the  primitive  elements  of  the  decomposition  are  flows  along  simple 
cycles  and  flows  along  simple  paths. 

Theorem  1.7.2  [26]  For  every  pseudoflow  /,  there  exists  a  collection  of  k  <  m  pseudoflows  gl, . . . 
such  that  for  every  i,  the  graph  induced  by  the  set  of  arcs  on  which  g ,  is  positive  is  either  a  simple 
cycle  or  a  simple  path  from  a  vertex  with  deficit  to  a  vertex  with  excess;  it  consists  of  arcs  (v,w)  with 
f(v,w)  >  0,  and  (1.10)  is  satisfied.  Such  a  decomposition  can  be  found  in  O(nm)  time. 

The  five  primitive  elements  of  the  decomposition  for  generalized  pseudoflows  g  are  defined  as 
follows.  The  set  of  arcs  on  which  an  element  of  the  decomposition  of  g  is  positive  is  a  subset  of  the 
arcs  on  which  g  is  positive.  Let  T(g)  denote  the  set  of  vertices  with  excesses  and  let  S{g)  denote 
the  set  of  vertices  with  deficits.  The  five  types  of  primitive  elements  are  classified  according  to  the 
graph  induced  by  the  set  of  arcs  with  positive  flow.  Type  I:  A  path  from  a  .crtex  in  S{g)  to  a 
vertex  in  T(g).  Such  a  flow  creates  a  deficit  and  an  excess  at  the  two  ends  of  the  path.  Type  II:  A 
flow-generating  cycle  and  a  path  leading  from  this  cycle  to  a  vertex  in  T(g).  Such  a  flow  creates 
excess  at  the  end  of  the  path.  (If  the  path  ends  at  the  source,  then  this  corresponds  to  a  GAP.) 
Type  III :  A  flow-absorbing  cycle  and  a  path  from  a  vertex  in  S(g)  to  this  cycle.  Such  a  flow  creates 
a  deficit  at  the  end  of  the  path.  Type  IV:  A  cycle  with  unit  gain.  Such  a  flow  does  not  create  any 
excesses  or  deficits.  Type  V:  A  pair  of  cycles  connected  by  a  path,  where  one  of  the  cycles  generates 
flow  and  the  other  one  absorbs  it.  Such  a  flow  does  not  create  any  excesses  or  deficits. 

Theorem  1.7.3  [42,  52]  For  a  generalized  pseudoflow  g  there  exist  k  <  m  primitive  pseudoflows 
gi,-..,gk  such  that,  for  each  i,  g<  is  of  one  of  the  five  types  described  above,  and  (110)  is  satisfied. 
Such  a  decomposition  can  be  found  in  0(nm )  time. 


Remark:  In  all  three  of  these  theorems  the  time  to  find  the  claimed  decomposition  can  be  reduced 
to  0(m  log  n)  by  using  the  dynamic  tree  data  structure  that  will  be  discussed  in  Section  2.o. 


Chapter  2 


The  Maximum  Flow  Problem 

2.1  Introduction 


The  maximum  flow  problem  has  been  studied  for  over  thirty  years.  The  classical  methods  for 
solving  this  problem  are  the  Ford- Fulkerson  augmenting  path  method  [25.  26],  the  closely  related 
blocking  flow  method  of  Dinic  [20],  and  an  appropriate  variant  of  the  network  simplex  method  (see 
e.g.  [60]).  A  push-relabel  method  for  solving  maximum  flow  problems  has  been  recently  proposed 
by  Goldberg  [39]  and  fully  developed  by  Goldberg  and  Tarjan  [44,  48].  Recently,  parallel  and 
distributed  algorithms  for  the  maximum  flow  problem  have  been  studied  as  well. 

Many  maximum  flow  algorithm  se  scaling  techniques  and  data  structures  to  achieve  better 
running  time  bounds.  The  scaling  techniques  used  by  maximum  flow  algorithms  are  capacity  scaling , 
introduced  by  Edmonds  and  Karp  [21]  in  the  context  of  the  minimum-cost  circulation  problem, 
and  the  closely  related  excess  scaling  of  Ahuja  and  Orlin  [3].  The  dynamic  tree  data  structure  of 
Sleator  and  Tarjan  [93,  94]  makes  it  possible  to  speed  up  many  maximum  flow  algorithms. 

The  technical  part  of  this  survey  deals  with  the  push-relabel  method,  which  gives  better  se¬ 
quential  and  parallel  complexity  bounds  than  previously  known  methods,  and  seems  to  outperform 
them  in  practice.  In  addition  to  describing  the  basic  method,  we  show  how  to  use  excess  scaling 
and  dynamic  trees  to  obtain  faster  implementations  of  the  method.  In  this  section  we  discuss  the 
previous  approaches  and  their  relationship  to  the  push-relabel  method. 

The  augmenting  path  algorithm  of  Ford  and  Fulkerson  is  based  on  the  fact  that  a  flow  is 
maximum  if  and  only  if  there  is  no  augmenting  path.  The  algorithm  repeatedly  finds  an  augmenting 
path  and  augments  along  it,  until  no  augmenting  path  exists.  Although  this  simple  generic  method 
need  not  terminate  if  the  network  capacities  are  reals,  and  it  can  run  in  exponential  time  if  the 
capacities  are  integers  represented  in  binary  [26],  it  can  be  made  .-fficient  by  restricting  the  choice 
of  augmenting  paths.  Edmonds  and  Karp  [21]  have  shown  that  two  strategies  for  selecting  the 
next  augmenting  path  give  polynomial  time  bounds.  The  first  such  strategy  is  to  select,  at  each 
iteration,  a  shortest  augmenting  path,  where  the  length  of  a  path  is  the  number  of  arcs  on  the  path. 
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The  second  strategy  is  to  select,  at  each  iteration,  a  fattest  path,  where  the  fatness  of  a  path  is  the 
minimum  of  the  residual  capacities  of  the  arcs  on  the  path. 

Independently,  Dinic  [20]  suggested  finding  augmenting  paths  in  phases,  handling  all  augmenting 
paths  of  a  given  shortest  length  in  one  phase  by  finding  a  blocking  flow  in  a  layered  network.  The 
shortest  augmenting  path  length  increases  from  phase  to  phase,  so  that  the  blocking  flow  method 
terminates  in  n  —  1  phases,  Dinic’s  discovery  motivated  a  number  of  algorithms  for  the  blocking 
flow  problem.  Dinic  proposed  an  O(nm)  blocking  flow  algorithm.  Soon  thereafter,  Karzanov  [68] 
proposed  an  0(n 2)  algorithm,  which  achieves  the  best  known  bound  for  dense  graphs.  Karzanov 
also  introduced  the  concept  of  a  preflow;  this  concept  is  used  by  many  maximum  flow  algorithms. 
Simpler  algorithms  achieving  an  0(n2)  bound  arc  described  in  [74,  97].  A  sequence  of  algorithms 
achieving  better  and  better  running  times  on  non-dense  graphs  has  been  proposed  [17,  32,  35. 
36,  46,  91,  94].  The  algorithm  of  [32]  uses  capacity  scaling;  the  algorithms  of  [94]  and  [46]  use 
dynamic  trees.  The  fastest  currently-known  sequential  algorithm  for  the  blocking  flow  problem, 
due  to  Goldberg  and  Tarjan,  runs  in  0(m  log(n2/m))  time  [46]. 

The  push-relabel  method  [39,  40,  44,  46]  replaces  the  layered  network  of  the  blocking  flow 
method  by  a  distance  labeling  of  the  vertices,  which  gives  an  estimate  of  the  distance  to  the  sink 

„he  lesiuual  graph.  The  algorithm  maintains  a  preflow  and  a  distance  labeling,  and  uses  two 
simple  operations,  pushing  and  relabeling ,  to  update  the  preflow  and  the  labeling,  repeating  them 
until  a  maximum  flow  is  found.  The  pushing  operation  is  implicit  in  Karzanov’s  algorithm.  The 
construction  of  a  layered  network  at  each  iteration  of  Dinic’s  method  can  be  viewed  as  a  sequence 
of  relabeling  operations.  Unlike  the  blocking  flow  method,  the  push-relabel  method  solves  the 
maximum  flow  problem  directly,  without  reducing  the  problem  to  a  sequence  of  blocking  flow 
subproblems.  As  a  result,  this  method  is  more  general  and  flexible  than  the  blocking  flow  method, 
and  leads  to  improved  sequential  and  parallel  algorithms.  Ahuja  and  Orlin  [3]  suggested  one  way 
of  recasting  algorithms  based  on  Dinic’s  approach  into  the  push-relabel  framework. 

Goldberg’s  FIFO 1  algorithm  [39]  runs  in  0{n3)  time.  Goldberg  and  Tarjan  [44,  48]  showed 
how  to  use  the  dynamic  tree  data  structure  to  improve  the  running  time  of  this  algorithm  to 
0(nm  log(n2/m)).  They  also  gave  a  largest-label  variant  of  the  algorithm,  for  which  they  obtained 
the  same  time  bounds  -  0{n3)  without  dynamic  trees  and  O(nmlog(n2/m))  with  such  trees. 
Cheriyan  and  Maheshwari  [16]  showed  that  (without  dynamic  trees)  the  FIFO  algorithm  runs 
in  0(n3)  time,  whereas  the  largest-label  algorithm  runs  in  0(n2y/m)  time.  Ahuja  and  Orlin  [3] 
described  an  0(nm  +  n2  log  f/)-time  version  of  the  push-relabel  method  based  on  excess  scaling. 
Ahuja,  Orlin,  and  Tarjan  [4]  gave  a  modification  of  the  Ahuja-Orlin  algorithm  that  runs  in  0(nm  + 
n2 v/log  U )  time  without  the  use  of  dynamic  trees  and  in  0(nm  log(^\/l0g  U  +  2))  time  with  them. 

The  primal  network  simplex  method  (see  e.g.  [60])  is  another  classical  method  for  solving  the 
maximum  flow  problem.  Cunningham  [18]  gave  a  simple  pivoting  rule  that  avoids  cycling  duo  to 

’FIFO  is  an  abbreviation  of  first-in,  first-out. 
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degeneracy.  Of  course,  if  a  simplex  algorithm  does  not  cycle,  it  must  terminate  in  exponential  time. 
Until  recently,  no  polynomial  time  pivoting  rule  was  known  for  the  primal  network  simplex  method. 
GoMfarb  and  Hao  [51]  have  given  such  a  rule.  The  resulting  algorithm  does  0(nm)  pivots,  which 
correspond  to  0(nm)  flow  augmentations;  it  runs  in  0(n2m)  time.  Interestingly,  the  blocking  flow 
method  also  makes  0(nm )  flow  augmentations.  Unlike  the  blocking  flow  method,  the  Goldfarb- 
Hao  algorithm  does  not  necessarily  augment  along  shortest  augmenting  paths.  In  their  analysis, 
Goldfarb  and  Hao  use  a  variant  of  distance  labeling  and  a  variant  of  the  relabeling  operation 
mentioned  above.  Dynamic  trees  can  be  used  to  obtain  an  0(nm  log  n)-time  implementation  of 
their  algorithm  [41]. 

The  best  currently  known  sequential  bounds  for  the  maximal  flow  problem  are  0{nm\og(n2 /m)) 
and  0(um  log(^  v^°g  U  +  2)).  Note  that,  although  the  running  times  of  the  known  algorithms  come 
very  close  to  an  0(mn)  bound,  the  existence  of  a  maximum  flow  algorithm  that  meets  this  bound 
remains  open. 

With  the  increasing  inheres*  in  parallel  computing,  parallel  and  distributed  algorithms  for  the 
maximum  flow  problem  have  received  a  great  deal  of  attention.  The  first  parallel  algorithm  for  the 
problem,  due  to  Shiloach  and  Vishkin  [92],  runs  in  0(n2logn)  time  on  an  n-processor  PRAM  [27] 
and  uses  0{n)  memory  per  processor.  In  a  synchronous  distributed  model  of  computation,  this  al¬ 
gorithm  runs  in  0(n2)  time  using  0(n3)  messages  and  0(n)  memory  per  processor.  The  algorithm 
of  Goldberg  [39,  40,  48]  uses  less  memory  than  that  of  Shiloach  and  Vishkin:  O(m)  memory  for 
the  PRAM  implementation  and  0(A)  memory  per  processor  for  the  distributed  implementation 
(where  A  is  the  processor  degree  in  the  network).  The  time,  processor,  and  message  bounds  of 
this  algorithm  are  the  same  as  those  of  the  Shiloach-Vishkin  algorithm.  Ahuja  and  Orlin  [3]  devel¬ 
oped  a  PRAM  implementation  of  their  excess-scaling  algorithm.  The  resource  bounds  are  [m/n] 
processors,  O(m)  memory,  and  0(n2  log  n  log  U)  time.  Chcriyan  and  Maheshwari  [16]  proposed 
a  synchronous  distributed  implementation  of  the  largest-label  algorithm  that  runs  in  0(n2)  time 
using  O(A)  memory  per  processor  and  0(n\/m)  messages. 

For  a  long  time,  the  primal  network  simplex  method  was  the  method  of  choice  in  practice.  A 
study  of  Goldfarb  and  Grigoriadis  [50]  suggested  that  the  algorithm  of  Dinic  [20]  performs  better 
than  the  network  simplex  method  and  better  than  the  later  algorithms  based  on  the  blocking  flow 
method.  Recent  studies  of  Ahuja  and  Orlin  (personal  communication)  and  Grigoriadis  (personal 
communication)  show  superiority  of  various  versions  of  the  push-relabel  method  to  Dinic's  algo¬ 
rithm.  An  experimental  study  of  Goldberg  [40]  shows  that  a  substantial  speedup  can  be  achieved 
by  implementing  the  FIFO  algorithm  on  a  highly  parallel  computer. 

More  efficient  algorithms  have  been  developed  for  the  special  case  of  planar  networks.  Ford  and 
Fulkerson  [25]  have  observed  that  the  the  maximum  flow  problem  on  a  planar  network  is  related 
to  a  shortest  path  problem  on  the  planar  dual  of  the  network.  The  algorithms  in  [8.  25,  98,  55,  56, 


Figure  2.1:  The  generic  maximum  flow  algorithm. 


59,  63,  75,  87]  make  clever  use  of  this  observation. 

2.2  A  Generic  Algorithm 

In  this  section  we  describe  the  generic  push-relabel  algorithm  [40,  44,  48].  First,  however,  we  need 
the  following  definition.  For  a  given  preflow  /,  a  distance  labeling  is  a  function  d  from  the  vertices 
to  the  non-negative  integers  such  that  d(t)  =  0,  d(s)  =  n,  and  d(v)  <  d(w)  +  1  for  all  residual  arcs 
(v,w).  The  intuition  behind  this  definition  is  as  follows.  Define  a  distance  graph  G‘j  as  follows. 
Add  an  arc  (s,t)  to  Gj.  Define  the  length  of  all  residual  arcs  to  be  equal  to  one  and  the  length  of 
the  arc  (s,<)  to  be  n.  Then  d  is  a  “locally  consistent”  estimate  on  the  distance  to  the  sink  in  the 
distance  graph.  (In  fact,  it  is  easy  to  show  that  d  is  a  lower  bound  on  the  distance  to  the  sink.) 
We  denote  by  do-(v,w)  the  distance  from  vertex  v  to  vertex  w  in  the  distance  graph.  The  generic 

algorithm  maintains  a  preflow  /  and  a  distance  labeling  d  for  /,  and  updates  /  and  d  using  pxtsh 
and  relabel  operations.  To  describe  these  operations,  we  need  the  following  definitions.  We  say  that 
a  vertex  v  is  active  if  v  £  {s,t}  and  e/(v)  >  0.  Note  that  a  preflow  /  is  a  flow  if  and  only  if  there 
are  no  active  vertices.  An  arc  (v,w)  is  admissible  if  {v,w)  €  Ej  and  d(v)  =  d(w)  +  1. 

The  algorithm  begins  with  the  preflow  /  that  is  equal  to  the  arc  capacity  on  each  arc  leaving 
the  source  and  zero  on  all  arcs  not  incident  to  the  source,  and  with  some  initial  labeling  d.  The 
algorithm  then  repetitively  performs,  in  any  order,  the  update  operations ,  push  and  relabel ,  described 
in  Figure  2.2.  When  there  are  no  active  vertices,  the  algorithm  terminates.  A  summary  of  the 
algorithm  appears  in  Figure  2.1. 
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push(v,  tv). 

Applicability:  v  is  active  and  ( v,w )  is  admissible. 

Action:  send  6  €  (0,  min(ey  (t>),  uj  (v,  tc))]  units  of  flow  from  v  to  w. 

relabel(v). 

Applicability:  either  s  or  t  is  reachable  from  v  in  Gj  and  Vtc  €  V'  uj(v,w)  =  0  or  d{w )  >  d(v). 
Action:  replace  d(v)  by  min(„,U0€  Ef{d(w)}  +  1. 


Figure  2.2:  The  update  operations.  The  pushing  operation  updates  the  preflow,  and  the  relabeling 
operation  updates  the  distance  labeling.  Except  for  the  excess  scaling  algorithm,  all  algorithms  discussed 
in  this  section  push  the  maximum  possible  amount  6  when  doing  a  push. 

The  update  operations  modify  the  preflow  /  and  the  labeling  d.  A  push  from  v  to  u •  increases 
f(v,w)  and  e j(w)  by  up  to  f  =  min{e/(v),  u j( i\  u')}.  and  decreases  f(u\  v)  and  fj(v)  by  the  same 
amount.  The  push  is  saturating  if  uj(i\w)  =  0  after  the  push  and  nonsaturating  otherwise.  A 
relabeling  of  v  sets  the  label  of  v  equal  to  the  largest  value  allowed  by  the  valid  labeling  constraints. 

There  is  one  part  of  the  algorithm  we  have  not  yet  specified:  the  choice  of  an  initial  labeling 
d.  The  simplest  choice  is  d(s)  -  n  and  d(  v)  =  0  for  r  €  V-  —  {s}.  A  more  accurate  choice  (indeed, 
the  most  accurate  possible  choice)  is  d(v)  -  dG-{v,t)  for  r  G  V\  where  /  is  the  initial  preflow.  The 

latter  labeling  can  be  computed  in  0(n/)  time  using  backwards  breadth-first  searches  from  the  sink 
and  from  the  source  in  the  residual  graph.  The  resource  bounds  we  shall  derive  for  the  algorithm  are 
correct  for  any  valid  initial  labeling.  To  simplify  the  proofs,  we  assume  that  the  algorithm  starts 
with  the  simple  labeling.  In  practice,  it  is  preferable  to  start  with  the  most  accurate  values  of 
the  distance  labels,  and  to  update  the  distance  labels  periodically  by  using  backward  breadth-first 
search. 

Remark:  By  giving  priority  to  relabeling  operations,  it  is  possible  to  maintain  the  following  invari¬ 
ant:  Just  before  a  push,  d  gives  the  exact  distance  to  t  in  the  distance  graph.  Furthermore,  it  is 
possible  to  implement  the  relabeling  operations  so  that  the  total  work  done  to  maintain  the  distance 
labels  is  O(nm)  (see  e.g.  [48]).  Since  the  running  time  bounds  derived  in  this  section  are  fi (nm). 
one  can  assume  that  the  relabeling  is  done  in  this  way.  In  practice,  however,  maintaining  exact 
distances  is  expensive;  a  better  solution  is  to  maintain  a  valid  distance  labeling  and  periodically 
update  it  to  the  exact  labeling. 

Next  we  turn  our  attention  to  the  correctness  and  termination  of  the  algorithm.  Our  proof  of 
correctness  is  based  on  Theorem  1.2.1.  The  following  lemma  is  important  in  the  analysis  of  the 
algorithm. 

Lemma  2.2.1  If  /  is  a  preflow  and  v  is  a  vertex  with  positive  excess,  then  the  source  s  is  reachable 
from  r  in  the  residual  graph  Gj. 


Using  this  lemma  and  induction  on  the  number  of  update  operations,  it  can  be  shown  that 
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one  of  the  two  update  operations  must  be  applicable  to  an  active  vertex,  and  that  the  operations 
maintain  a  valid  distance  labeling  and  preflow. 

Theorem  2.2.2  [48]  Suppose  that  the  algorithm  terminates.  Then  the  preflow  /  is  a  maximum  flow. 

Proof :  When  the  algorithm  terminates,  ail  vertices  in  V  —  {s,t}  must  have  zero  excess,  because 
there  are  no  active  vertices.  Therefore  /  must  be  a  flow.  We  show  that  if  /  is  a  preflow  and  d  is 
a  valid  labeling  for  /,  then  the  sink  t  is  not  reachable  from  the  source  s  in  the  residual  graph  G /. 
Then  Theorem  1.2.1  implies  that  the  algorithm  terminates  with  a  maximum  flow. 

Assume  by  way  of  contradiction  that  there  is  an  augmenting  path  s  =  i’0.  t’i , . . . ,  i’/  =  t.  Then 
I  <  n  and  (ty,  r1+1 )  E  E;  for  0  <  i  <  l.  Since  d  is  a  valid  labeling,  we  have  d{  jy)  <  d{  r,  +  l )  -f  1  for 
0  <  i  <  l.  Therefore,  we  have  d(i  )  <  d(t)  +  l  <  n,  since  d(t)  —  0.  which  contradicts  d(s)  =  n.  | 

The  key  to  the  running  time  analysis  of  the  algorithm  is  the  following  lemma,  which  shows  that 
distance  labels  cannot  increase  too  much. 

Lemma  2.2.3  At  any  time  during  ihe  execution  of  the  algorithm,  for  any  vertex  v  6  V',  d(v)  <  2n  -  1. 

Proof:  The  lemma  is  trivial  for  v  =  s  and  v  =  t.  Suppose  r  G  V  -  {.«./},  Since  the  algorithm 
changes  vertex  labels  by  only  by  means  of  the  relabeling  operation,  it  is  enough  to  prove  the  lemma 
for  a  vertex  r  such  that  s  or  t  is  reachable  from  r  in  Gj.  Thus  there  is  a  simple  path  from  r  to  .« 
or  t  in  G j.  Let  r  =  t’o,  t'i , . . . ,  be  such  a  path.  The  length  /  of  the  path  is  at  most  n  -  1.  Since  d 
is  a  valid  labeling  and  ( r, .  rl+i )  E  Ej.  we  have  d{v, )  <  d(  r1  +  1 )  +  1 .  Therefore,  since  d{  rt )  is  either 
n  or  0.  we  have  d(  v)  =  d(  fo)  <  d{  r,  )  -f  /  <  n  +  (n  -  1 )  =  2n  -  1.  | 

Lemma  2.2.3  limits  the  number  of  relabeling  operations,  and  allows  us  to  amortize  the  work 
done  by  the  algorithm  over  increases  in  vertex  labels.  The  next  two  lemmas  bound  the  number  of 
relabelings  and  the  number  of  saturating  pushes. 

Lemma  2.2.4  The  number  of  relabeling  operations  is  at  most  2 n  —  1  per  vertex  and  at  most  (2 n  - 
l)(n  -  2)  <  2 n2  overall. 

Lemma  2.2.5  The  number  of  saturating  pushes  is  at  most  nm. 

Proof:  For  an  arc  (t\tr)  G  E ,  consider  the  saturating  pushes  from  r  to  u\  After  one  such  push. 
u/(v.u')  =  0,  and  another  such  push  cannot  occur  until  d(u-)  increases  by  at  least  2,  a  push  from  i r 
to  r  occurs,  and  d{ v)  increases  by  at  least  2.  If  we  charge  each  saturating  push  from  v  to  w  except 
the  first  to  the  preceding  label  increase  of  r,  we  obtain  an  upper  bound  of  n  on  the  number  of  such 
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The  most  interesting  part  of  the  analysis  is  obtaining  a  bound  on  the  number  of  nonsaturating 
pushes.  For  this  we  use  amortized  analysis  and  in  particular  the  potential  function  technique  (see 
e.g.  [98]). 

Lemma  2.2.6  The  number  of  nonsaturating  pushing  operations  is  at  most  2n2m. 

Proof:  We  define  the  potential  $  of  the  current  preflow  /  and  labeling  d  by  the  formula  $  = 
is  active}  d(v).  We  have  0  <  $  <  2 n2  by  Lemma  2.2.3.  Each  nonsaturating  push,  say  from  a 
vertex  v  to  a  vertex  w,  decreases  $  by  at  least  one,  since  d{w)  —  d(v)  —  1  and  the  push  makes  v 
inactive.  It  follows  that  the  total  number  of  nonsaturating  pushes  over  the  entire  algorithm  is  at 
most  the  sum  of  the  increases  in  <f>  during  the  course  of  the  algorithm,  since  $  =  0  both  at  the 
beginning  and  at  the  end  of  the  computation.  Increasing  the  label  of  a  vertex  v  by  an  amount 
k  increases  by  k.  The  total  of  such  increases  over  the  algorithm  is  at  most  2n2.  A  saturating 
push  can  increase  <h  by  at  most  2n  -  2.  The  total  of  such  increases  over  the  entire  algorithm  is  at 
most  (2 n  —  2)nm.  Summing  gives  a  bound  of  at  most  2 n2  +  (2 n  -  2 )nm  <  2 n2m  on  the  number  of 
nonsaturating  pushes.  | 


Theorem  2.2.7  [48]  The  generic  algorithm  terminates  after  0{n2m)  update  operations. 

Proof:  Immediate  from  Lemmas  2.2.4,  2.2.5.  and  2.2.6.  | 

The  running  time  of  the  generic  algorithm  depends  upon  the  order  in  which  update  operations 
are  applied  and  on  implementation  details.  In  the  next  sections  we  explore  these  issues.  First  we  give 
a  simple  implementation  of  the  generic  algorithm  in  which  the  time  required  for  the  nonsaturating 
pushes  dominates  the  overall  running  time.  Sections  2.3  and  2.4  specify  orders  of  t he  update 
operations  that  decrease  the  number  of  nonsaturating  pushes  and  permit  O ( rr 3 )  and  0(mn  t 
:  2  log  f/)-time  implementations.  Section  2.5  explores  an  orthogonal  approach.  It  shows  how  to  use 
sophisticated  data  structures  to  reduce  the  time  per  nonsaturating  push  rather  than  the  number 
of  such  pushes. 

2.3  Efficient  Implementation 

Our  first  step  toward  an  efficient  implementation  is  a  way  of  combining  the  update  operations 
locally.  We  need  some  data  structures  to  represent  the  network  and  the  preflow.  We  call  an 
unordered  pair  {r,u*}  such  that  (r,tr)  £  £  an  edge  of  G.  We  associate  the  three  values  n(v.u-). 
u(u’,r),  and  /(r,u’)(  =  —  /(tr,r))  with  each  edge  {t\u-}.  Each  vertex  r  has  a  list  of  the  incident 
edges  {r,  tr},  in  fixed  but  arbitrary  order.  Thus  each  edge  {t’.tc}  appears  in  exactly  two  lists,  the 
one  for  v  and  the  one  for  u'.  Each  vertex  v  has  a  current  edge  { t\  ir).  which  is  the  current  candidate 
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discharge(v). 

Applicability:  v  is  active. 

Action:  let  be  the  current  edge  of  r; 

itme-io-relabel  *—  false; 
repeat 

if  (n,  ic)  is  admissible  then  pvsh(v,w) 
else 

if  {v,  u>)  is  not  the  last  edge  on  the  edge  list  of  t;  then 

replace  {v,  w}  as  the  current  edge  of  n  by  the  next  edge  on  the  list 
else  begin 

make  the  first  edge  on  the  edge  li;  of  the  current  edge; 
ttme-to-relabel  <—  true; 
end; 

until  ej( v)  =  0  or  time-to-relabel; 
if  iime-to-relabel  then  relabel{v); 

Figure  2.3:  The  discharge  operation. 

for  a  pushing  operation  from  v.  Initially,  the  current  edge  of  v  is  the  first  edge  on  the  edge  list  of 
r.  The  main  loop  of  the  implementation  consists  of  repeating  the  discharge  operation  described  in 
Figure  2.3  until  there  are  no  active  vertices.  (We  shall  discuss  the  maintenance  of  active  vertices 
later.)  The  discharge  operation  is  applicable  to  an  active  vertex  v.  This  operation  iteratively 
attempts  to  push  the  excess  at  v  through  the  current  edge  of  v  if  a  pushing  operation  is 

applicable  to  this  edge.  If  not,  the  operation  replaces  {u,u>}  as  the  current  edge  of  v  by  the  next 
edge  on  the  edge  list  of  v;  or,  if  {u,u>}  is  the  last  edge  on  this  list,  it  makes  the  first  edge  on  the 
list  the  current  one  and  relabels  v.  The  operation  stops  when  the  excess  at  v  is  reduced  to  zero  or 
v  is  relabeled. 

The  following  lemma  shows  that  discharge  does  relabeling  correctly;  the  proof  of  the  lemma  is 
straightforward. 

Lemma  2.3.1  The  discharge  operation  does  a  relabefing  only  when  the  relabeling  operation  is  appli¬ 
cable. 

Lemma  2.3.2  The  version  of  the  generic  push/relabel  algorithm  based  on  discharging  runs  in  O(nm) 
time  plus  the  total  time  needed  to  do  the  nonsaturating  pushes  and  to  maintain  the  set  of  active  vertices. 

Any  representation  of  the  set  of  active  vertices  that  allows  insertion,  deletion,  and  access  to  some 
active  vertex  in  0(1)  time  results  in  an  0(n2m)  running  time  for  the  discharge- based  algorithm, 
by  Lemmas  2.2.6  ana  2.3.2.  (Pushes  can  be  implemented  in  0(1)  time  per  push.) 

By  processing  active  vertices  in  a  more  restricted  order,  we  obtain  improved  performance.  Two 
natural  orders  were  suggested  in  [44,  48].  One,  the  FIFO  algorithm ,  is  to  maintain  the  set  of  active 
vertices  as  a  queue,  always  selecting  for  discharging  the  front  vertex  on  the  queue  .  1  adding  newly 
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procedure  process-verier, 

remove  a  vertex  v  from  Bj; 
old-label «—  d( v); 
discharge^ ); 

add  each  vertex  w  made  active  by  the  discharge  to  Bj(Wy, 

if  d(v)  ^  old-label  then  begin 
b  —  d(t>); 
add  v  to  Bj; 

end 

else  if  Bt  =  0  then  6  <—  6  —  1; 
end. 


Figure  2.4:  The  process-vertex  procedure. 


active  vertices  to  the  rear  of  the  queue.  The  other,  the  largest-label  algorithm,  is  to  always  select 
for  discharging  a  vertex  with  the  largest  label.  The  FIFO  algorithm  runs  in  0(n3)  time  [41.  4K] 
and  the  largest-label  algorithm  runs  in  0(n2y/m)  time  [16].  We  shall  derive  an  0(n3)  time  bound 
for  both  algorithms,  after  first  describing  in  a  little  more  detail  how  to  implement  largest-label 
selection. 

The  implementation  maintains  an  array  of  sets  B,,  0  <  i  <  2n  -  1,  and  an  index  b  into  the 
array.  Set  B,  consists  of  all  active  vertices  with  label  i,  represented  as  a  doubly-linked  list,  so 
that  insertion  and  deletion  take  0(1)  time.  The  index  b  is  Die  largest  label  of  an  active  vertex. 
During  the  initialization,  when  the  arcs  going  out  of  the  source  are  saturated,  the  resulting  active 
vertices  are  placed  in  B0,  and  b  is  set  to  0.  At  each  iteration,  the  algorithm  removes  a  vertex 
from  B&,  processes  it  using  the  discharge  operation,  and  updates  b.  The  algorithm  terminates 
when  b  becomes  negative,  i.e.,  when  there  are  no  active  vertices.  This  processing  of  vertices,  which 
implements  the  while  loop  of  the  generic  algorithm,  is  described  in  Figure  2.4. 

To  understand  why  the  process-vertex  procedure  correctly  maintains  b.  note  that  dischargc(v) 
either  relabels  v  or  gets  rid  of  all  excess  at  r,  but  not  both.  In  the  former  case,  r  is  the  active 
vertex  with  the  largest  distance  label,  so  6  must  be  increased  to  d(v).  In  the  latter  case,  the  excess 
at  v  has  been  moved  to  vertices  with  distance  labels  of  b  -  1,  so  if  B/,  is  empty,  then  6  must  be 
decreased  by  one.  The  total  time  spent  updating  b  during  the  course  of  the  algorithm  is  0(n2). 

The  bottleneck  in  both  the  FIFO  method  and  the  largest-label  method  is  the  number  of  non¬ 
saturating  pushes.  We  shall  obtain  an  0(n3)  bound  on  the  number  of  such  pushes  by  dividing  the 
computation  into  phases,  defined  somewhat  differently  for  each  method.  For  the  FIFO  method, 
phase  1  consists  of  the  discharge  operations  applied  to  the  vertices  added  to  the  queue  by  the 
initialization  of  /;  phase  i  +  1,  for  i  <  1,  consists  of  the  discharge  operations  applied  to  the  vertices 
added  to  the  queue  during  phase  i.  For  the  largest-label  method,  a  phase  consists  of  a  maximal 
interval  of  time  during  which  b  remains  constant. 


Lemma  2.3.3  The  number  of  phases  during  the  running  of  either  the  FIFO  or  the  largest-label  algo- 
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Proof:  Define  the  potential  $  of  the  current  /  and  d  by  $  =  max|v|„  ;8  active}d(i’)>  with  the  max¬ 
imum  taken  to  be  zero  if  there  are  no  active  vertices.  (In  the  case  of  the  largest-label  algorithm, 
$  =  b  except  on  termination.)  There  can  be  only  2n2  phases  that  do  one  or  more  relabelings.  A 
phase  that  does  no  relabeling  decreases  $  by  at  least  one.  The  initial  and  final  values  of  $  are  zero. 
Thus  the  number  of  phases  that  do  no  relabeling  is  at  most  the  sum  of  the  increases  in  4>  during 
the  computation.  The  only  increases  in  $  axe  due  to  label  increases;  an  increase  of  a  label  bv  k 
can  cause  $  to  increase  by  up  to  k.  Thus  the  sum  of  the  increases  in  $  over  the  computation  is  at 
most  2n2,  and  so  is  the  number  of  phases  that  do  no  relabeling.  | 


Theorem  2.3.4  [48]  Both  the  FIFO  and  the  largest-label  algorithm  run  in  0(n3)  time. 

Proof:  For  either  algorithm,  there  is  at  most  one  nonsaturating  push  per  vertex  per  phase.  Thus 
by  Lemma  2.3.3  the  total  number  of  nonsaturating  pushes  is  0(n3),  as  is  the  running  time  by 
Lemma  2.3.2.  | 

Cheriyan  and  Maheshwari  [16],  by  means  of  an  elegant  balancing  argument,  were  able  to  improve 
the  bound  on  the  number  of  nonsaturating  pushes  in  the  largest-label  algorithm  to  0{  n2v/m),  giving 
the  following  result: 

Theorem  2.3.5  [16]  The  largest-label  algorithm  runs  in  0{n2y/m)  time. 

2.4  Excess  Scaling 

A  different  approach  to  active  vertex  selection  leads  to  running  time  bounds  dependent  on  the  size 
U  of  the  largest  capacity  as  well  as  on  the  graph  size.  This  approach,  excess  scaling ,  was  introduced 
by  Ahuja  and  Orlin  [3]  and  developed  further  by  Ahuja,  Orlin,  and  Tarjan  [4].  We  shall  describe 
in  detail  a  slight  revision  of  the  original  excess-scaling  algorithm,  which  has  a  running  time  of 
0(nm  +  n2log(/). 

For  the  termination  of  the  excess-scaling  method,  all  arc  capacities  must  be  integral;  hence 
we  assume  throughout  this  section  that  this  is  the  case.  The  method  preserves  integrality  of  the 
flow  throughout  the  computation.  It  depends  on  a  parameter  A  that  is  an  upper  bound  on  the 
maximum  excess  of  an  active  vertex.  Initially  A  =  2^og^ .  The  algorithm  proceeds  in  phases;  after 
each  phase,  A  is  halved.  When  A  <  1,  all  active  vertex  excesses  must  be  zero,  and  the  algorithm 
terminates.  Thus  the  algorithm  terminates  after  at  most  log 2U  +  1  phases.  To  maintain  the 
invariant  that  no  active  vertex  excess  exceeds  A,  the  algorithm  does  not  always  push  the  maximum 
possible  amount  when  doing  a  pushing  operation.  Specifically,  when  pushing  from  a  vertex  r  to  a 
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vertex  w,  the  algorithm  moves  an  amount  of  flow  6  given  by  6  =  min{e/(i;),  uj(v,w ),  A  -  ej(w)} 
if  w  /  t,6  =  min{ey(v),  tiy(r,tc)}  if  w  =  t.  That  is,  6  is  the  maximum  amount  that  can  be  pushed 
while  maintaining  the  invariant. 

The  algorithm  consists  of  initialization  of  the  preflow,  the  distance  labels,  and  A,  followed  by 
a  sequence  of  process-vertex  operations  of  the  kind  described  in  Section  2.3.  Vertex  selection  for 
process-vertex  operations  is  done  by  the  large  excess,  smallest  label  rule:  process  an  active  vertex  v 
with  e/(v)  >  A/2;  among  all  candidate  vertices,  choose  one  of  smallest  label.  When  every  active 
vertex  v  has  e/(v)  <  A/2,  A  is  halved  and  a  new  phase  begins;  when  there  are  no  active  vertices, 
the  algorithm  stops. 

Since  the  excess- scaling  algorithm  is  a  version  of  the  generic  process-vertex- based  algorithm 
described  in  Section  2.3,  Lemma  2.3.2  applies.  The  following  lemma  bounds  the  number  of  nonsat¬ 
urating  pushes: 


Lemma  2.4.1  The  number  of  nonsaturating  pushes  during  the  excess-scaling  algorithm  is  0(n2  log  U). 

Proof:  We  define  a  potential  $  by  «f>  =  is  active}  ef(v)d(v)/ A.  Since  e/(v) /A  <  1  for  any 

active  vertex,  0  <  $  <  2n2.  Every  pushing  operation  decreases  <J>.  Consider  a  nonsaturating 
push  from  a  vertex  v  to  a  vertex  tv.  The  large  excess,  smallest  label  rule  guarantees  that  before 
the  push  ej( v)  >  A/2  and  either  ej(w)  <  A/2  or  w  =  t.  Thus  the  push  moves  at  least  A/2 
units  of  flow,  and  hence  decreases  $  by  at  least  1/2.  The  initial  and  final  values  of  are  zero,  so 
the  total  number  of  nonsaturating  pushes  is  at  most  twice  the  sum  of  the  increases  in  <f>  over  the 
course  of  the  algorithm.  Increasing  the  label  of  a  vertex  by  k  can  increase  $  by  at  most  k.  Thus 
relabelings  account  for  a  total  increase  in  $  of  at  most  2 n2.  A  change  in  phase  also  increases  4>,  by 
a  factor  of  two,  or  at  most  n2 .  Hence  the  sum  of  increases  in  $  is  at  most  2 n2  +  n 2  (log2  V  +  1). 
and  therefore  the  number  of  nonsaturating  pushes  is  at  most  An2  +  2n2(log2  U  +  1).  | 

The  large  excess,  smallest  label  rule  can  be  implemented  by  storing  the  active  vertices  with 
excess  exceeding  A/2  in  an  array  of  sets,  as  was  previously  described  for  the  largest-label  rule.  In 
the  former  case,  the  index  6  indicates  the  nonempty  set  of  smallest  index.  Since  b  decreases  only 
when  a  push  occurs,  and  then  by  at  most  1,  the  total  time  spent  updating  b  is  O  ( nm  +  n2  logf’). 
From  Lemmas  2.3.2  and  2.4.1  we  obtain  the  following  result: 


Theorem  2.4.2  [3]  The  excess- sea  ling  algorithm  runs  in  0  (nm  +  n2  log  U)  time. 


A  more  complicated  version  of  excess  scaling,  devised  by  Ahuja,  Orlin,  and  Tarjan  [4],  has 
a  running  time  of  O  (nm  +  n2\/log  U).  This  algorithm  uses  a  hybrid  vertex  selection  rule  that 
combines  a  stack-based  mechanism  with  the  “wave”  approach  of  Tarjan  [97]. 
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Figure  2.5:  The  dynamic  tree  operations. 

2.5  Use  of  Dynamic  Trees 

The  two  previous  sections  discussed  ways  of  reducing  the  number  of  nonsaturating  pushes  by 
restricting  the  order  of  the  update  operations.  An  orthogonal  approach  is  to  Teduce  the  time  per 
nonsaturating  push  rather  than  the  number  of  such  pushes.  The  idea  is  to  perform  a  succession 
of  pushes  along  a  single  path  in  one  operation,  using  a  sophisticated  data  structure  to  make  this 
possible.  Observe  that  immediately  after  a  nonsaturating  push  along  an  arc  (v,w),  {v,w)  is  still 
admissible,  and  we  know  its  residual  capacity.  The  dynamic  tree  data  structure  of  Sleator  and 
Tarjan  [93,  94]  provides  an  efficient  way  to  maintain  information  about  such  arcs.  We  shall  describe 
a  dynamic  tree  version  of  the  generic  algorithm  that  has  a  running  time  of  O  ( nm  log  n). 

The  dynamic  tree  data  structure  allows  the  maintenance  of  a  collection  of  vertex-disjoint  rooted 
trees,  each  arc  of  which  has  an  associated  real  value.  We  regard  a  tree  arc  as  directed  toward  the 
root,  t.e.,  from  child  to  parent.  We  denote  the  parent  of  a  vertex  v  by  parent(v),  and  adopt  the 
convention  that  every  vertex  is  both  an  ancestor  and  a  descendant  of  itself.  The  data  structure 
supports  the  seven  operations  described  in  Figure  2.5.  A  sequence  of  /  tree  operations  on  trees  of 
maximum  size  (number  of  vertices)  k  takes  O(flogfc)  time. 

In  the  dynamic  tree  algorithm,  the  arcs  of  the  dynamic  trees  axe  a  subset  of  the  admissible  arcs 
and  every  active  vertex  is  a  tree  root.  The  value  of  an  arc  (v,w)  is  its  residual  capacity.  Initially 
each  vertex  is  made  into  a  one-vertex  dynamic  tree  using  a  make-tree  operation. 

The  heart  of  the  dynamic  tree  implementation  is  the  tree-push  operation  described  in  Figure  2.6. 
This  operation  is  applicable  to  an  admissible  arc  (v,w)  such  that  v  is  active.  The  operation  adds 
(v,w)  to  the  forest  of  dynamic  trees,  pushes  as  much  flow  as  possible  along  the  tree  from  v  to  the 
root  of  the  tree  containing  w,  and  deletes  from  the  forest  each  arc  that  is  saturated  by  this  flow 


make-tree(v):  Make  vertex  v  into  a  one-vertex  dynamic  tree.  Vertex  v  must  be  in  no  other  tree. 

find-root(v):  Find  and  return  the  root  of  the  tree  containing  vertex  v. 

find-value(v):  Find  and  return  the  value  of  the  tree  arc  connecting  v  to  its  parent.  If  v  is  a  tree 
root,  return  infinity. 

find-min(v):  Find  and  return  the  ancestor  w  of  v  such  that  the  tree  arc  connecting  w  to  its  parent 
has  minimum  value  along  the  path  from  v  to  find-root(v).  In  case  of  a  tie,  choose 
the  vertex  w  closest  to  the  tree  root.  If  t;  is  a  tree  root,  return  v. 


change-value(v,  x): 
link{v,  w,  x): 

cut(v): 


Add  real  number  x  to  the  value  of  every  arc  dong  the  path  from  t>  to  find-rooi(v). 

Combine  the  trees  containing  v  and  w  by  making  w  the  parent  of  v  and  giving  the 
new  tree  arc  joining  v  and  ui  the  value  x.  This  operation  does  nothing  if  v  and  w 
are  in  the  same  tree  or  if  v  is  not  a  tree  root. 

Break  the  tree  containing  v  into  two  trees  by  deleting  the  arc  from  t’  to  its  parent. 
This  operation  does  nothing  if  v  is  a  tree  root. 
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Tree-push  (e,tc) 

Applicability:  v  is  active  and  (e.te)  is  admissible. 

Action:  link  (t>,  w,  u/(v,  w)); 

parent  (v)  «—  w; 

6  <—  min  {ej  (u),  find-value(find-min(v))‘, 
change-value  (t 

while  v  ^  find-rool(v)  and  find-value(find-min(v))  =  0  do  begin 
z  *—  find-min(v); 
cui  (c); 

f(z,  parenl(z))  *—  u(z,  parent (z)); 
f  (parent  (z),z)  * - u(z,  parent(z)); 

end. 


Figure  2.6:  The  tree-push  operation. 


change. 

The  dynamic  tree  algorithm  is  just  the  generic  algorithm  with  the  push  operation  replaced  by 
the  tree-push  operation,  with  the  initialization  modified  to  make  each  vertex  into  a  dynamic  tree, 
and  with  a  postprocessing  step  that  extracts  the  correct  flow  value  for  each  arc  remaining  in  a 
dynamic  tree.  (This  postprocessing  takes  one  find-value  operation  per  vertex.) 

It  is  straightforward  to  show  that  the  dynamic  tree  implementation  is  correct,  by  observing  that 
it  maintains  the  invariant  that  between  tree-pushing  operations  every  active  vertex  is  a  dynamic 
tree  root  and  every  dynamic  tree  arc  is  admissible.  The  following  lemma  bounds  the  number  of 
tree-push  operations: 

Lemma  2.5.1  The  number  of  the  tree-push  operations  done  by  the  dynamic  tree  implementation  is 
O(nm). 

Proof:  Each  tree-push  operation  either  saturates  an  arc  (thus  doing  a  saturating  push)  or  decreases 
the  number  of  active  vertices  by  one.  The  number  of  active  vertices  can  increase,  by  at  most  one, 
as  a  result  of  a  tree-push  operation  that  saturates  an  arc.  By  Lemma  2.2.5  there  are  at  most  nm 
saturating  pushes.  An  upper  bound  of  2 nm  +  n  on  the  number  of  tree-push  operations  follows, 
since  initially  there  are  at  most  n  active  vertices.  | 

If  we  implement  the  dynamic  tree  algorithm  using  the  discharge  operation  of  Section  2.3  with 
push  replaced  by  tree-push,  we  obtain  the  following  result  (assuming  active  vertex  selection  takes 
0(1)  time): 

Theorem  2.5.2  The  discharge- based  implementation  of  the  dynamic  tree  algorithm  runs  in  0(nm  log  n) 
time. 


Proof:  Each  tree-pushing  operation  does  0(1)  dynamic  tree  operations  plus  0(1)  per  arc  saturated. 
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The  theorem  follows  from  Lemma  2.2.5,  Lemma  2.3.2,  and  Lemma  2.5.1,  since  the  maximum  size 
of  any  dynamic  tree  is  0(n).  | 

The  dynamic  tree  data  structure  can  be  used  in  combination  with  the  FIFO,  largest-label, 
or  excess-scaling  method.  The  resulting  time  bounds  are  0(nm\og(n2 / m))  for  the  FIFO  and 
largest-label  methods  [48]  and  0(nm log(^  i/log  U)  +  2))  for  the  fastest  version  of  the  excess- 
scaling  method  [4].  In  each  case,  the  dynamic  tree  method  must  be  modified  so  that  the  trees  are 
not  too  large,  and  the  analysis  needed  to  obtain  the  claimed  bound  is  rather  complicated. 


Chapter  3 

The  Minimum- Cost  Circulation 
Problem:  Cost- Scaling 

3.1  Introduction 


Most  polynomial-time  algorithms  for  the  minimum-cost  circulation  problem  use  the  idea  of  scaling. 
This  idea  was  introduced  by  Edmonds  and  Karp  [21],  who  used  it  to  develop  the  first  polynomial¬ 
time  algorithm  for  the  problem.  Scaling  algorithms  work  by  obtaining  a  sequence  of  feasible  or 
almost-feasible  solutions  that  are  closer  and  closer  to  the  optimum.  The  Edmonds-Karp  algorithm 
scales  capacities.  Rock  [89]  was  the  first  to  propose  an  algorithm  that  scales  costs.  Later,  Bland  and 
Jensen  [13]  proposed  a  somewhat  different  cost-scaling  algorithm,  which  is  closer  to  the  generalized 
cost  scaling  method  discussed  in  this  chapter. 

The  cost-scaling  approach  works  by  solving  a  sequence  of  problems,  Fq,  Fj,  . . . .  Ft.  on  the  net¬ 
work  with  original  capacities  but  approximate  costs.  The  cost  function  c,  for  F,  is  obtained  by 
taking  the  i  most  significant  bits  of  the  original  cost  function  c.  The  first  problem  F0  has  zero 
costs,  and  therefore  the  zero  circulation  is  optimal.  An  optimal  solution  to  problem  F,_i  can  be 
used  to  obtain  an  optimal  solution  to  problem  F;  in  at  most  n  maximum  flow  computations  [13,  89]. 
Note  that  for  k  =  fJog2C],  c*  =  c.  Thus  the  algorithm  terminates  in  O(nlogC)  maximum  flow- 
computations. 

Goldberg  and  Tarjan  [40,  45,  47]  proposed  a  generalized  cost-scaling  approach.  The  idea  of 
thi6  method  (which  is  described  in  detail  below)  is  to  view  the  maximum  amount  of  violation  of 
the  complementary  slackness  conditions  as  an  error  parameter,  and  to  improve  this  parameter  by 
a  factor  of  tw’o  at  each  iteration.  Initially  the  error  is  at  most  C,  and  if  the  error  is  less  than 
1/n,  then  the  current  solution  is  optimal.  Thus  the  generalized  cost-scaling  method  terminates 
in  0(\ag(nC))  iterations.  The  computation  done  at  each  iteration  is  similar  to  a  maximum  flow 
computation.  The  traditional  cost-scaling  method  of  Rock  also  improves  the  error  from  iteration  to 
iteration,  bi  t  it  does  so  indirectly,  by  increasing  the  precision  of  the  current  costs  and  solving  the 
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resulting  problem  exactly.  Keeping  the  original  costs,  as  does  the  generalized  cost-scaling  approach, 
makes  it  possible  to  reduce  the  number  of  iterations  required  and  to  obtain  strongly-polynomial 
running  time  bounds.  Chapter  4  discusses  a  strongly-polynomial  version  of  the  generalized  cost¬ 
scaling  method.  For  further  discussion  of  generalized  versus  traditional  cost-scaling,  see  [46]. 

Time  bounds  for  the  cost-scaling  algorithms  mentioned  above  are  as  follows.  The  algorithms  of 
Rock  [89]  and  Bland  and  Jensen  [13]  run  in  0(n  log(C)M(n,  m,  U))  time,  where  M(n,m,U )  is  the 
time  required  to  compute  a  maximum  flow  on  a  network  with  n  vertices,  m  arcs,  and  maximum  arc 
capacity  U.  As  we  have  seen  in  Chapter  2,  M  =  0(nm  log min{n2/m,  -^N/log  U  4-  2}).  The  fastest 
known  implementation  of  the  generalized  cost-scaling  method  runs  in  0(nmlog(n2/m)log(nC)) 
time  [45].  It  is  possible  to  combine  cost  scaling  with  capacity  scaling.  The  first  algorithm  that 
combines  the  two  scaling  techniques  is  due  to  Gabow  and  Tarjan  [34].  A  different  algorithm  was 
proposed  by  Ahuja  et  al.  [1].  The  latter  algorithm  runs  in  0(nm  log  log  l/iog(nC))  time,  which 
makes  it  the  fastest  known  algorithm  for  the  problem  under  the  similarity  assumption. 


3.2  Approximate  Optimality 

A  key  notion  is  that  of  approximate  optimality ,  obtained  by  relaxing  the  complementary  slackness 
constraints  in  Theorem  1.3.3.  For  a  constant  e  >  0,  a  pseudoflow  /  is  said  to  be  (-optimal  u'ith 
respect  to  a  price  function  p  if,  for  every  arc  ( v ,  w),  we  have 

f(v,w)  <  u(v,w)  =>  cp(v,w)  >  -6  (e-optimality  constraint).  (3.1) 

A  pseudoflow  /  is  (-optimal  if  /  is  e-optimal  with  respect  to  some  price  function  p. 

An  important  property  of  e-optimality  is  that  if  the  arc  costs  are  integers  and  e  is  small  enough, 
any  e-optimal  circulation  is  minimum-cost.  The  following  theorem,  of  Bertsekas  [11]  captures  this 
fact. 


Theorem  3.2.1  [11]  If  ail  costs  are  integers  and  e  <  l/n,  then  an  e-optimal  circulation  /  is  minimum- 
cost. 


The  e-optimality  constraints  were  first  published  by  Tardos  [95]  in  a  paper  describing  the  first 
strongly  polynomial  algorithm  for  the  minimum-cost  circulation  problem.  Bertsekas  [11]  proposed 
a  pseudopolynomial  algorithm  based  upon  Theorem  3.2.1;  his  algorithm  makes  use  of  a  fixed 
e  <  l/n.  Goldberg  and  Tarjan  [40,  46,  47]  devised  a  successive  approximation  scheme  that  produces 
a  sequence  of  circulations  that  axe  e-optimal  for  smaller  and  smaller  values  of  c;  when  e  is  small 
enough,  the  scheme  terminates  with  an  optimal  circulation.  We  discuss  this  scheme  below. 


3.3.  THE  GENERALIZED  COST-SCALING  FRAMEWORK 
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procedure  min-cost(V,  E,u,c)- 

[initialization] 

e*-C; 

Vr,  p(v)  *-  0; 

V(v,  u>)  €  E,  f(v,  w)  <—  0; 

[loop] 

while  (  >  1/n  do 

—  refine(c,f,p); 

return(/); 

end. 


Figure  3.1:  The  generalized  cost-scaling  method. 

3.3  The  Generalized  Cost-Scaling  Framework 

Throughout  the  rest  of  this  chapter,  we  assume  that  all  arc  costs  are  integral.  We  give  here  a  high- 
level  description  of  the  generalized  cost-scaling  method  (see  Figure  3.1).  The  algorithm  maintains 
an  error  parameter  e,  a  circulation  /  and  a  price  function  p,  such  that  /  is  c-optimal  with  respect  to 
p.  The  algorithm  starts  with  e  =  C  (or  alternatively  c  =  2flog2cl),  with  p(v)  =  0  for  all  v  6  V,  and 
with  the  zero  circulation.  Any  circulation  is  C-optimal.  The  main  loop  of  the  algorithm  repeatedly 
reduces  the  error  parameter  e.  When  c  <  1/n,  the  current  circulation  is  minimum-cost,  and  the 
algorithm  terminates. 

The  task  of  the  subroutine  refine  is  to  reduce  the  error  in  the  optimality  of  the  current  circu¬ 
lation.  The  input  to  refine  is  an  error  parameter  c,  a  circulation  /,  and  a  price  function  p  such 
that  /  is  c-optimal  with  respect  to  p.  The  output  from  refine  is  a  reduced  error  parameter  c,  a 
new  circulation  /,  and  a  new  price  function  p  such  that  /  is  c-optimal  with  respect  to  p.  The 
implementations  of  refine  described  in  this  survey  reduce  the  error  parameter  c  by  a  factor  of  two. 

The  correctness  of  the  algorithm  is  immediate  from  Theorem  3.2.1,  assuming  that  refine  is 
correct.  The  number  of  iterations  of  refine  is  0(log(nC)).  This  gives  us  the  following  theorem: 

Theorem  3.3.1  [47]  A  minimum-cost  circulation  can  be  computed  in  the  time  required  for 
0(log(nC))  iterations  of  refine,  if  refine  reduces  c  by  a  factor  of  at  least  two. 

3.4  A  Generic  Refinement  Algorithm 

In  this  section  we  describe  an  implementation  of  refine  that  is  a  common  generalization  of  the 
generic  maximum  flow  algorithm  of  Chapter  2.2  and  the  auction  algorithm  for  the  assignment 
problem  [9]  (first  published  in  [10]).  We  call  this  the  generic  implementation.  This  implementation, 
proposed  by  Goldberg  and  Tarjan  [47],  is  essentially  the  same  as  the  main  loop  of  the  minimum- 
cost  circulation  algorithm  of  Bertsekas  [11],  which  is  also  a  common  generalization  of  the  maximum 
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procedure  refine(e,  f,p); 

[initialization] 
c  —  e/2; 

V(t),ui)  £  E  do  if  cp(v,w)  <  0  then  begin  f(v,w)  «—  u(v,w),  f(w,v)  < - u(v,w);  end; 

[loop] 

while  3  an  update  operation  that  applies  do 
select  such  an  operation  and  apply  it; 
return(e,  /,  p); 

end. 


Figure  3.2:  The  generic  refine  subroutine. 


push(v  ,w). 

Applicability:  v  is  active  and  (v,w)  is  admissible. 

Action:  send  6  =  min(ey  (t>),  uj  (v,  xv))  units  of  flow  from  t>  to  w. 

rtlabel(v). 

Applicability:  v  is  active  and  Vic  £  V  {u/(v,  w )  =  0  or  cp(v,  w)  >  0 ). 
Action:  replace  p(v)  by  max(t)  u,)££;  {p(tu)  -  c(t>,  te)  -  e). 


Figure  3.3:  The  update  operations  for  the  generic  refinement  algorithm.  Compare  with  Figure  2.2. 

flow  and  assignment  algorithms.  The  ideas  behind  the  auction  algorithm  can  be  used  to  give  an 
alternative  interpretation  to  the  results  of  [40,  47]  in  terms  of  relaxation  methods;  see  [12]. 

As  we  have  mentioned  in  Section  3.3,  the  effect  of  refine  is  to  reduce  f  by  a  factor  of  two  while 
maintaining  the  e-optimality  of  the  current  flow  /  with  respect  to  the  current  price  function  p. 
The  generic  refine  subroutine  is  described  on  Figure  3.2.  It  begins  by  halving  e  and  saturating 
every  arc  with  negative  reduced  cost.  This  converts  the  circulation  /  into  an  e-optimal  pseudoflow 
(indeed,  into  a  0-optimal  pseudoflow).  Then  the  subroutine  converts  the  e-optimal  pseuuofiow  into 
an  e-optimal  circulation  by  applying  a  sequence  of  the  update  operations  push  and  relabel ,  each  of 
which  preserves  e-optimality. 

The  inner  loop  of  the  generic  algorithm  consists  of  repeatedly  applying  the  two  update  op¬ 
erations,  described  in  Figure  3.3,  in  any  order,  until  no  such  operation  applies.  To  define  these 
operations,  we  need  to  redefine  admissible  arcs  in  the  context  of  the  minimum-cost  circulation 
problem.  Given  a  pseudoflow  /  and  a  price  function  p,  we  say  that  an  arc  ( v,w )  is  admissible  if 
(v,  w)  is  a  residual  arc  with  negative  reduced  cost. 

A  push  operation  applies  to  an  admissible  arc  ( v,w )  such  that  vertex  v  is  active.  It  consists  of 
pushing  6  =  min{e/(u),  uj(v,  w)}  units  of  flow  from  v  to  u\  thereby  decreasing  e/(v)  and  /(«•,  r) 
by  6  and  increasing  e/(te)  and  f(v,w)  by  6.  The  push  is  saturating  if  uj(v,w)  =  0  after  the  push 
and  nonsaturating  otherwise. 
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A  relabel  operation  applies  to  an  active  vertex  v  that  has  no  exiting  admissible  arcs.  It  con¬ 
sists  of  decreasing  p(v)  to  the  smallest  value  allowed  by  the  e-optimalitv  constraints,  namely 
max(„>u,)€E/ {-c(w,  to)  +  p(w)  -  e}. 

If  an  e-optimal  pseudoflow  /  is  not  a  circulation,  then  either  a  pushing  or  a  relabeling  operation 
is  applicable.  It  is  easy  to  show  that  any  pushing  operation  preserves  e-optimality.  The  next  lemma 
gives  two  important  properties  of  the  relabeling  operation. 

Lemma  3.4.1  Suppose  /  is  an  e-optimal  pseudoflow  with  respect  to  a  price  function  p  and  a  vertex  v 
is  relabeled.  Then  the  price  of  t>  decreases  by  at  least  e  and  the  pseudoflow  /  is  e-optimal  with  respect 
to  the  new  price  function  p' . 

Proof:  Before  the  relabeling,  cp(v,tv)  >  0  for  all  (u,u>)  G  £/,  i.e.,  p(v)  >  p(ie)  -  cp(v,w )  for  all 
{v ,«-)  G  Ej.  Thus  p'(v)  =  max(ViVl)€£/  {p(w)  -  c(v,w )  -  e)  <  p(v)  -  e. 

To  verify  e-optimality,  observe  that  the  only  residual  arcs  whose  reduced  costs  are  affected  by 
the  relabeling  are  those  of  the  form  (v,  w)  or  (w,  v).  Any  arc  of  the  form  (u\v)  has  its  reduced  cost 
increased  by  the  relabeling,  preserving  its  e-optimality  constraint.  Consider  a  residual  arc  (v.ir). 
By  the  definition  of  p',  p'(t’)  >  p(ir)  -  c(t>, u>)  —  e.  Thus  cp>(v,w)  =  c(v.u')  +  p'(v)  -  p(ir)  >  — e, 
which  means  that  (v,  tr)  satisfies  its  e-optimality  istraint.  | 

Since  the  update  operations  preserve  e-optimality,  and  since  some  update  operation  applies  if  / 
is  not  a  circulation,  it  follows  that  if  refine  terminates  and  returns  (c,/,p),  then  /  is  a  circulation 
which  is  e-optimal  with  respect  to  p.  Thus  refine  is  correct. 

Next  we  analyze  the  number  of  update  operations  that  can  take  place  during  an  execution  of 
refine.  We  begin  with  a  definition.  The  admissible  graph  is  the  graph  Ga  —  (T,  Ea)  such  that 
Ea  is  the  set  of  admissible  arcs.  As  refine  executes,  the  admissible  graph  changes.  An  important 
invariant  is  that  the  admissible  graph  remains  acyclic. 

Lemma  3.4.2  Immediately  after  a  relabeling  is  applied  to  a  vertex  v,  no  admissible  arcs  enter  v. 

Proof  :  Let  (u,r)  be  a  residual  arc.  Before  the  relabeling,  cp(u,  r)  >  — <  by  c-optimality.  By 
Lemma  3.4.1,  the  relabeling  decreases  p( r),  and  hence  increases  cp(u.v ),  by  at  least  c.  Thus 
cp(u,v)  >  0  after  the  relabeling.  | 

Corollary  3.4.3  Throughout  the  running  of  refine,  the  admissible  graph  is  acyclic. 

Proof :  Initially  the  admissible  graph  contains  no  arcs  and  is  thus  acyclic.  Pushes  obviously  preserve 
acyclicity.  Lemma  3.4.2  implies  that  reiabelings  also  preserve  acyclicity.  | 


Next  we  derive  a  crucial  lemma,  which  generalizes  Lemma  2.2.1. 
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Lemma  3.4.4  Let  /  be  a  pseudoflow  and  /'  a  circulation.  For  any  vertex  v  with  ej(v)  >  0,  there  is 
a  vertex  w  with  ej(w)  <  0  and  a  sequence  of  distinct  vertices  v  —  t’o,  tq , . . . ,  v/_i,  t ’/  =  w  such  that 
(v,,t;,+i)  £  Ej  and  (v.+i,^)  £  Ej<  for  0  <  t  <  /. 

Proof:  Let  v  be  a  vertex  with  ej{v)  >  0.  Define  G+  =  (V,  £+),  where  E+  =  {(x,y)\f'(x.y)  > 
f{x,y)},  and  define  G-  =  (V,E-),  where  E-  =  {(x,y)\f(x,y)  >  f'(x,y)}.  Then  E+  C  Ej,  since 
(x,y)  £  E+  implies  f(x,y)  <  f(x,y)  <  u(x,y).  Similarly  E-  C  Ep.  Furthermore  ( x,y )  £  E+  if 
and  only  if  ( y,x )  £  E_  by  antisymmetry.  Thus  to  prove  the  lemma  it  suffices  to  show  the  existence 
in  G+  of  a  simple  path  v  =  t>0,  vi  . . . ,  vt  with  e/(t y)  <  0. 

Let  5  be  the  set  of  all  vertices  reachable  from  v  in  G+  and  let  S  =  V  —  S.  (Set  S  may  be 
empty.)  For  every  vertex  pair  ( x,y )  £  S  x  S,  /(^,t/)  >  f'(x,y),  for  otherwise  y  £  S.  We  have 

0  =  E(x,y)e(SxS)n£  /'(*,  y)  since  /'  ^  a  circulation 

<  E(x,y)e(SxS)n£  Z(x’  V )  holds  term-by-term 

=  £(x,j,)e(SxS)nE/(*,2/)  +  E(x,v)€{Sx5)n£/(x-2/)  by  antisymmetry 

=  E(ity)e(SxV)n£  bv  definition  of  S 

=  -ExeSe/(x)  by  antisymmetry. 

But  v  £  S.  Since  e/(t>)  >  0,  some  vertex  «’  6  5  must  have  ey(tr)  <  0.  | 

Using  Lemma  3.4.4  we  can  bound  the  amount  by  which  a  vertex  price  can  decrease  during  an 
invocation  of  refine. 

Lemma  3.4.5  The  price  of  any  vertex  r  decreases  by  at  most  3 ne  during  an  execution  of  refine. 

Proof  :  Let  f2(  and  p 2,  be  the  circulation  and  price  functions  on  entry  to  refine.  Suppose  a  relabeling 
causes  the  price  of  a  vertex  v  to  decrease.  Let  /  be  the  pseudoflow  and  p  the  price  function  just 
after  the  relabeling.  Then  Cf(v)  >  0.  Let  v  =  r0,  tq, . . . ,  iy  =  w  with  ej(u-)  <  0  be  the  vertex- 
sequence  satisfying  Lemma  3.4.4  for  /  and  f  =  f2(. 

The  c-optimality  of  /  implies 


l-i  i-\ 

~  l(  <  X>p(r„rt+1)  =  p{ v)  -  p(w)  +  ^c(r,,r,+ ,). 

•=o  »=o 

The  2c-optimaJity  of  f2(  implies 


(3.2) 


/-i  ;-i 

~  2U  ^  5ZCP2t(r'  +  l'r')  =  P2r(«’)  -  P2r(l')  +  £1  C(  l’«  +  l '  r<  )• 

i=0  1=0 


(3.3) 


3.5.  EFFICIENT  IMPLEMENTATION 


39 


But  c(«i, v,-+i )  =  —  5Zi=oc(l,*+i’w«)  ^v  cost  antisymmetry.  Furthermore,  p(u-)  =  P2<(u') 

since  during  refine,  the  initialization  step  is  the  only  one  that  makes  the  excess  of  some  vertices 
negative,  and  a  vertex  with  negative  excess  has  the  same  price  as  long  as  its  excess  remains  negative. 
Adding  inequalities  (3.2)  and  (3.3)  and  rearranging  terms  thus  gives 

p(v)  >  P2({v)  -  3/e  >  P2c(v)  -  3 ne.  | 

Now  we  count  update  operations.  The  following  lemmas  are  analogous  to  Lemmas  2.2.4,  2.2.5 
and  2.2.6. 

Lemma  3.4.6  The  number  of  relabelings  during  an  execution  of  refine  is  at  most  3n  per  vertex  and 
3n(n  —  1)  in  total. 

Lemma  3.4.7  The  number  of  saturating  pushes  during  an  execution  of  refine  is  at  most  3nm. 

Proof:  For  an  arc  (v,w),  consider  the  saturating  pushes  along  this  arc.  Before  the  first  such  push 
can  occur,  vertex  v  must  be  relabeled.  After  such  a  push  occurs,  r  must  be  relabeled  before  another 
such  push  can  occur.  But  v  is  relabeled  at  most  3n  times.  Summing  over  all  arcs  gives  the  desired 
bound.  | 


Lemma  3.4.8  The  number  of  nonsaturating  pushes  during  one  execution  of  refine  is  at  most  3 r>2(  m  -f 

n). 

Proof:  For  each  vertex  r,  let  4>(r)  be  the  number  of  vertices  reachable  from  v  in  the  current 
admissible  graph  Ga ■  Let  <t>  =  0  if  there  are  no  active  vertices,  and  let  4>  =  ]T{<J>(r)|r  is  active} 
otherwise.  Throughout  the  running  of  refine ,  4>  >  0.  InitiaJly  4>  <  n.  since  G  a  has  no  arcs. 

Consider  the  effect  on  4>  of  update  operations.  A  nonsaturating  push  decreases  4>  by  at  least  one. 
since  G a  is  always  acyclic  by  Corollary  3.4.3.  A  saturating  push  can  increase  $  by  at  most  n.  since  at 
most  one  inactive  vertex  becomes  active.  If  a  vertex  v  is  relabeled,  4>  also  can  increase  by  at  most  », 
since  4>(u’)  for  tr  /  v  can  only  decrease  bv  Lemma  3.4.2.  The  total  number  of  nonsaturating  pushes 
is  thus  bounded  by  the  initial  value  of  plus  the  total  increase  in  <E>  throughout  the  algorithm. 
i.e.,  by  n  +  3 n2(n  —  1)  -f  3 n2m  <  3 n2(m  +  n).  | 


3.5  Efficient  Implementation 

As  in  the  case  of  the  generic  maximum  flow  algorithm,  we  can  obtain  an  especially  efficient  version 
of  refine  by  choosing  the  order  of  the  update  operations  carefully. 

Local  ordering  of  the  basic  operations  is  achieved  using  the  data  structures  and  the  discharge 
operation  of  Section  2.3.  The  discharge  operation  is  the  same  as  the  one  in  Figure  2.3,  but  uses  the 
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procedure  first-active; 

let  I,  be  a  list  of  all  vertices; 
let  v  be  the  first  vertex  in  L; 
while  3  an  active  vertex  do  begin 
if  v  is  active  then  begin 
discharge(v); 

if  the  discharge  has  relabeled  v  then 
move  v  to  the  front  of  L; 

end; 

else  replace  v  by  the  vertex  after  t>  on  L, 

end; 

end. 


Figure  3.4:  The  first-active  method. 

minimum-cost  circulation  versions  of  the  update  operations  (see  Figure  3.3).  As  in  the  maximum 
flow  case,  it  is  easy  to  show  that  discharge  applies  the  relabeling  operation  correctly.  The  overall 
running  time  of  refine  is  0{nm)  plus  0(1)  per  nonsaturating  push  plus  the  time  needed  for  active 
vertex  selection. 

In  the  maximum  flow  case,  one  of  the  vertex  selection  methods  we  considered  was  the  largest- 
label  method.  In  the  minimum-cost  circulation  case,  a  good  method  is  first-active  [46],  which 
is  a  generalization  of  the  largest-label  method.  The  idea  of  the  method  is  to  process  vertices  in 
topological  order  with  respect  to  the  admissible  graph. 

The  first-active  method  maintains  a  list  L  of  all  the  vertices  of  G,  in  topological  order  with 
respect  to  the  current  admissible  graph  G'a,  t.e.,  if  (i\tr)  is  an  arc  of  Ga •  v  appears  before  tr  in  L. 
Initially  L  contains  the  vertices  of  G  in  any  order.  The  method  consists  of  repeating  the  following 
step  until  there  are  no  active  vertices:  Find  the  first  active  vertex  on  L,  say  v,  apply  a  discharge 
operation  to  v,  and  move  v  to  the  front  of  L  if  the  discharge  operation  has  relabeled  r. 

In  order  to  implement  this  method,  we  maintain  a  current  vertex  v  of  L.  which  is  the  next 
candidate  for  discharging.  Initially  c  is  the  first  vertex  on  L.  The  implementation,  described  in 
Figure  3.4,  repeats  the  following  step  until  there  are  no  active  vertices:  If  v  is  active,  apply  a 
discharge  operation  to  it,  and  if  this  operation  relabels  u,  move  v  to  the  front  of  L;  otherwise  (t.e., 
v  is  inactive),  replace  v  by  the  vertex  currently  after  it  on  L.  Because  the  reordering  of  L  maintains 
a  topological  ordering  with  respect  to  G^,  no  active  vertex  precedes  v  on  L.  This  implies  that  the 
implementation  is  correct. 

Define  a  phase  as  a  period  of  time  that  begins  with  r  equal  to  the  first  vertex  on  L  and  ends 
when  the  next  relabeling  is  performed  (or  when  the  algorithm  terminates). 

Lemma  3.5.1  The  first-active  method  terminates  after  0{n2)  phases. 


Proof:  Each  phase  except  the  last  ends  with  a  relabeling  operation.  | 
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Theorem  3.5.2  [47]  The  first-active  implementation  of  refine  runs  in  0(n3)  time,  giving  an 
0(n3log(nC))  bound  for  finding  a  minimum-cost  circulation. 

Proof :  Lemma  3.5.1  implies  that  there  are  0(n3)  nonsaturating  pushes  (one  per  vertex  per  phase) 
during  an  execution  of  refine.  The  time  spent  manipulating  L  is  0(n)  per  phase,  for  a  total  of 
0(n3).  All  other  operations  require  a  total  of  0(nm)  time.  | 

A  closely  related  strategy  for  selecting  the  next  vertex  to  process  is  the  wave  method  [40. 
46,  47],  which  gives  the  same  0(n3)  running  time  bound  for  refine.  (A  similar  pseudopolynomial 
algorithm,  without  the  use  of  scaling  and  missing  some  of  the  implementation  details,  was  developed 
independently  in  [11].)  The  only  difference  between  the  first-active  method  and  the  wave  method 
is  that  the  latter,  after  moving  a  vertex  v  to  the  front  of  L ,  replaces  v  by  the  vertex  after  the  old 
position  of  v;  if  v  is  the  last  vertex  on  L,  v  is  replaced  by  the  first  vertex  on  L. 

As  in  the  maximum  flow  case,  the  dynamic  tree  data  structure  can  be  used  to  obtain  faster 
implementations  of  refine.  A  dynamic  tree  implementation  of  the  generic  version  of  refine  analogous 
to  the  maximum  flow  algorithm  discussed  in  Section  2.5  runs  in  0(nm  logn)  time  [47],  A  dynamic 
tree  implementation  of  either  the  first-active  method  or  the  wave  method  runs  in  0(nm  log(n2/m)) 
time  [46].  In  the  latter  implementation,  a  second  data  structure  is  needed  to  maintain  the  list  L. 
The  details  are  somewhat  involved. 

3.6  Refinement  Using  Blocking  Flows 

An  alternative  way  to  implement  the  refine  subroutine  is  to  generalize  Dinic’s  approach  to  the 
maximum  flow  problem.  Goldberg  and  Tarjan  [46,  47]  showed  that  refinement  can  be  carried  out 
by  solving  a  sequence  of  0(n)  blocking  flow  problems  on  acyclic  networks  (i.c.,  on  networks  for 
which  the  residual  graph  of  the  zero  flow  is  acyclic);  this  extends  Dinic's  result,  which  reduces  a 
maximum  flow  problem  to  n  -  1  blocking  flow  problems  on  layered  networks.  In  this  section  we 
describe  the  Goldberg-Tarjan  algorithm.  At  the  end  of  this  section,  we  make  a  few  comments  about 
blocking  flow  algorithms. 

To  describe  the  blocking  flow  version  of  refine  we  need  some  standard  definitions.  Consider 
a  flow  network  (G,u,s,t).  A  flow  /  is  blocking  if  any  path  from  s  to  t  in  the  residual  graph  of 
zero  flow  contains  a  saturated  arc,  i.e.,  an  arc  (v,w)  such  that  uj(v,w)  =  0.  A  maximum  flow 
is  blocking,  but  not  conversely.  A  directed  graph  is  layered  if  its  vertices  can  be  assigned  integer 
layers  in  such  a  way  that  layerfv )  =  layer{u>)  +  1  for  every  arc  (t>,  w).  A  layered  graph  is  acyclic 
but  not  conversely. 

An  observation  that  is  crucial  to  this  section  is  as  follows.  Suppose  we  have  a  pseudoflow  / 
and  a  price  function  p  such  that  the  vertices  can  be  partitioned  into  two  sets,  5  and  5,  such  that 
no  admissible  arc  leads  from  a  vertex  in  S  to  a  vertex  in  5;  in  other  words,  for  every  residual  arc 
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procedure  rtfine(c,  f,p); 

[initialization] 

f  *- f/2; 

V(t i,w)  €  E  do  if  cp(v,w)  <  0  then  f(v,w)  *—  u(v,tv); 

[loop] 

while  /  is  not  a  circulation  do  begin 

5  «—  {t;  G  Vj3u  €  V  such  that  ej(u)  >  0  and  v  is  reachable  from  u  in  G^}: 

Ve  €  5,  p(v)  —  p(t>)  -  e; 

let  N  be  the  network  formed  from  Ga  by  adding  a  source  s,  a  sink  f, 
an  arc  (s,t>)  of  capacity  ej(v)  for  each  v  €  V  with  e/(v)  >  0,  and 
an  arc  (v,t)  of  capacity  —ej{v)  for  each  v  €  V  with  tj(v)  <  0; 
find  a  blocking  flow  b  on  N; 

V(v,  tv)  €  EA,  f( v,  w)  *—  f(v,  w)  +  b(v,w); 

end; 

return  (c,/,p); 

end. 

Figure  3.5:  The  blocking  refine  subroutine. 

(v,w)  6  Ej  such  that  v  6  5  and  tv  G  5.  we  have  cp(v,w )  >  0.  Define  p'  to  be  equal  to  p  on  S 
and  to  p  -  e  on  S.  It  is  easy  to  see  that  replacing  p  by  p'  does  not  create  any  new'  residual  arc 
with  reduced  cost  less  than  —  e.  The  blocking  flow  method  augments  by  a  blocking  flow  to  create 
a  partitioning  of  vertices  as  described  above,  and  modifies  the  price  function  by  replacing  p  by  p'. 

Figure  3.5  describes  an  implementation  of  refine  that  reduces  c  by  a  factor  of  two  by  computing 
0{n )  blocking  flows.  This  implementation  reduces  e  by  a  factor  of  two,  saturates  all  admissible 
arcs,  and  then  modifies  the  resulting  pseudoflow  (while  maintaining  e-optimality  with  respect  to 
the  current  price  function)  until  it  is  a  circulation.  To  modify  the  pseudoflow,  the  method  first 
partitions  the  vertices  of  G  into  two  sets  5  and  5,  such  that  5  contains  all  vertices  reachable  in 
the  current  admissible  graph  GA  from  vertices  of  positive  excess.  Vertices  in  5  have  their  prices 
decreased  by  c.  Next,  an  auxiliary  network  N  is  constructed  by  adding  to  G a  a  source  s,  a  sink  t , 
an  arc  (s,v)  of  capacity  ej{v)  for  each  vertex  v  with  e/(v)  >  0,  and  an  arc  (v,t)  of  capacity  — e/(r) 
for  each  vertex  with  e/(v)  <  0.  An  arc  (v,w)  €  EA  has  capacity  u/(v,w>)  in  N.  A  blocking  flow  b 
on  N  is  found.  Finally,  the  pseudoflow  /  is  replaced  by  the  pseudoflow  f'(v,w)  =  f(v,w)  +  b{  tv,  tv) 
for  ( v ,  tv)  €  E. 

The  correctness  of  the  blocking  flow  method  follows  from  the  next  lemma,  which  can  be  proved 
by  induction  on  the  number  of  iterations  of  the  method. 

Lemma  3.6.1  The  set  S  computed  in  the  inner  loop  contains  only  vertices  v  with  ej(  v)  >  0.  At  the 
beginning  of  an  iteration  of  the  loop,  /  is  an  c-optimal  pseudoflow  with  respect  to  the  price  function 
p.  Decreasing  the  prices  of  vertices  in  S  by  e  preserves  the  c-optimality  of  /.  The  admissible  graph 
remains  acyclic  throughout  the  algorithm. 
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The  bound  on  the  number  of  iterations  of  the  method  follows  from  Lemma  3.4.5  Bind  the  fact 
that  the  prices  of  the  vertices  with  deficit  remain  unchanged,  while  the  prices  of  the  vertices  with 
excess  decrease  by  c  during  every  iteration. 

Lemma  3.0.2  The  number  of  iterations  of  the  inner  loop  in  the  blocking  flow  implementation  of  refine 
is  at  most  3  n. 

Since  the  running  time  of  an  iteration  of  the  blocking  flow  method  is  dominated  by  the  time 
needed  for  a  blocking  flow  computation,  we  have  the  following  theorem. 

Theorem  3.6.3  [47]  The  blocking  flow  implementation  of  refine  runs  in  0(nB(n,m))  time,  giving  an 
0(nB(n,  m)log(nC))  bound  for  finding  a  minimum-cost  circulation,  where  B(n,m)  is  the  time  needed 
to  find  a  blocking  flow  in  an  acyclic  network  with  n  vertices  and  m  arcs. 

The  fastest  known  sequential  algorithm  for  finding  a  blocking  flow  on  an  acyclic  network  is  due 
to  Goldberg  and  Tarjan  [46]  and  runs  in  0(mlog(n2/m))  time.  Thus,  by  Theorem  3.6.3,  we  obtain 
an  0( nm log(n2/m)  log( nC))  time  bound  for  the  minimum-cost  circulation  problem.  This  is  the 
same  as  the  fastest  known  implementation  of  the  generic  refinement  method. 

There  is  a  crucial  difference  between  Dinic’s  maximum  flow  algorithm  and  the  blocking  flow- 
version  of  refine.  Whereas  the  former  finds  blocking  flows  in  layered  networks,  the  latter  must  find 
blocking  flows  in  acyclic  networks,  an  apparently  harder  task.  Although  for  sequential  computation 
the  acyclic  case  seems  to  be  no  harder  than  the  layered  case  (the  best  know-n  time  bound  is 
0(mlog(n2/m)  for  both),  this  is  not  true  for  sequential  computation.  The  Shiloach-Vishkin  PRAM 
blocking  flow  algorithm  [92]  for  layered  networks  runs  in  0(n2  logn)  time  using  0(n2)  memory  and 
n  processors.  The  fastest  known  PRAM  algorithm  for  acyclic  networks,  due  to  Goldberg  and 
Tarjan  [49],  runs  in  the  same  0(n2  \ogn)  time  bound  but  uses  0{nm)  memory  and  m  processors. 
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Chapter  4 

Strongly  Polynomial  Algorithms 
Based  on  Cost-Scaling 


4.1  Introduction 

The  question  of  whether  the  minimum-cost  circulation  problem  has  a  strongly  polynomial  algorithm 
was  posed  in  1972  by  Edmonds  and  Karp  [21]  and  resolved  only  in  1984  by  Tardos  [95].  Her  result 
led  to  the  discovery  of  a  number  of  strongly  polynomial  algorithms  for  the  problem  [30,  37,  82].  In 
this  chapter  we  discuss  several  strongly  polynomial  algorithms  based  on  cost  scaling;  in  the  next, 
we  explore  capacity-scaling  algorithms,  including  strongly  polynomial  ones.  Of  the  known  strongly 
polynomial  algorithms,  the  asymptotically  fastest  is  the  capacity-scaling  algorithm  of  Orlin  [82]. 

We  begin  by  describing  a  modification  of  the  generalized  cost-scaling  method  that  makes  it 
strongly  polynomial  [46].  Then  we  describe  the  minimum-mean  cycle-canceling  algorithm  of  Gold¬ 
berg  and  Tarjan  [45].  This  simple  algorithm  is  a  specialization  of  Klein’s  cycle-canceling  algorithm 
[70];  it  does  not  use  scaling,  but  its  analysis  relies  on  ideas  related  to  cost  scaling. 


4.2  Fitting  Price  Functions  and  Tight  Error  Parameters 


In  order  to  obtain  strongly  polynomial  bounds  on  the  generalized  cost-scaling  method,  we  need  to 
take  a  closer  look  at  the  notion  of  e-optimality  defined  in  Section  3.2.  The  definition  of  e-optimality 
motivates  the  following  two  problems; 


1.  Given  a  pseudoflow  /  and  a  constant  e  >  0,  find  a  price  function  p  such  that  /  is  e-optimal 
with  respect  to  p,  or  show  that  there  is  no  such  price  function  (i.e.,  that  /  is  not  e-optimal). 

2.  Given  a  pseudoflow  /,  find  the  the  smallest  e  >  0  such  that  /  is  e-optimal.  For  this  e.  we  say 
that  /  is  c-tight. 
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The  problem  of  finding  an  optima]  price  function  given  an  optimal  circulation  is  the  special  case 
of  Problem  1  with  e  =  0.  We  shall  see  that  the  first  problem  can  be  reduced  to  a  shortest  path 
problem,  and  that  the  second  problem  requires  the  computation  of  a  cycle  of  minimum  average  arc 
cost. 

To  address  these  problems.  we  need  some  results  about  shortest  paths  and  shortest  pat1  .re.' 
(see  e.g.  [96]).  Let  G  be  a  directed  graph  with  a  distinguished  source  vertex  s  from  which  every  vertex 
is  reachable  and  a  cost  c(v,w)  on  every  arc  (v,w).  For  a  spanning  tree  T  rooted  at  s,  the  tree  cost 
function  d  :  V  — ►  R  is  defined  recursively  as  follows:  d(s)  =  0,  d(v)  =  d(parent(v))  +  c(parent(v),v) 
for  v  6  V  —  {s},  where  parent(v)  is  the  parent  of  u  in  T.  A  spanning  tree  T  rooted  at  s  is  a  shortest 
path  tree  if  and  only  if,  for  every  vertex  t>,  the  path  from  s  to  v  in  T  is  a  minimum-cost  path  from 
s  to  v  in  G,  t.e.,  d(v)  is  the  cost  of  a  minimum-cost  path  from  s  to  v. 

Lemma  4.2.1  (see  e.g.  [96])  Graph  G  contains  a  shortest  path  tree  if  and  only  if  G  does  not  contain  a 
negative-cost  cycle.  A  spanning  tree  T  rooted  at  s  is  a  shortest  path  tree  if  and  only  if  c{v,w)  +  d(v)  > 
d{w)  for  every  arc  ( v,w )  in  G. 

Consider  Problem  1:  given  a  pseudoflow  /  and  a  nonnegative  c,  find  a  price  function  p  with 
respect  to  which  /  is  e-optimal,  or  show  that  /  is  not  e-optimal.  Define  a  new  cost  function 
c(e>  :  E  —  R  by  c<«>(v  ,  w)  =  c(v,w)  4-  e.  Extend  the  residual  graph  Gf  by  adding  a  single  vertex  s 
and  arcs  from  it  to  all  other  vertices  to  form  an  auxiliary  graph  Gaux  =  (Vaux,£aux)  =  (V  U  {s}, 
Ej  U  ({s}  x  V)).  Extend  c^  to  GOUx  by  defining  c^(s,  v)  =  0  for  every  arc  (s,i>),  where  v  e  V. 
Note  that  every  vertex  is  reachable  from  s  in  Gaux. 

Theorem  4.2.2  Pseudoflow  /  is  c-optimal  if  and  only  if  Gour  (or  equivalently  G j)  contains  no  cycle 
of  negative  c^-cost.  If  T  is  any  shortest  path  tree  of  Gaux  (rooted  at  s)  with  respect  to  the  arc  cost 
function  c'(\  and  d  is  the  associated  tree  cost  function,  then  /  is  c-optimal  with  respect  to  the  price 
function  p  defined  by  p(v)  =  d(v)  for  all  v  6  V. 

Proof:  Suppose  /  is  c-optimal.  Any  cycle  in  Gaux  is  a  cycle  in  Gj ,  since  vertex  s  has  no  incoming 
arcs.  Let  T  be  a  cycle  of  length  l  in  Gaux.  Then  c(T)  >  — fc,  which  implies  c^(T)  =  c(r)  +  U  >  0. 
Therefore  Gaux  contains  no  cycle  of  negative  c(f)-cost. 

Suppose  Gaux  contains  no  cycle  of  negative  c^-cost.  Then,  by  Lemma  4.2.1,  Gaux  has  some 
shortest  path  tree  rooted  at  s.  Let  T  be  any  such  tree  and  let  d  be  the  tree  cost  function.  By 
Lemma  4.2.1,  c^(v,w)  +  d(v)  >  d(w)  for  all  (v,w)  €  Ej,  which  is  equivalent  to  c(v,w)  +  d(v)-d(u') 
>  -c  for  all  (v,w)  €  Ej.  But  these  are  the  e-optimality  constraints  for  the  price  function  p  =  d. 
Thus  /  is  c-optimal  with  respect  to  p.  | 
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Using  Theorem  4.2.2,  we  can  solve  Problem  1  by  constructing  the  auxiliary  graph  Gaxix  and 
finding  either  a  shortest  path  tree  or  a  negative-cost  cycle.  Constructing  Gaux  takes  0(m )  time. 
Finding  a  shortest  path  tree  or  a  negative-cost  cycle  takes  0(nm )  time  using  the  Bellman-Ford 
shortest  path  algorithm  (see  e.g.  [96]). 

Let  us  turn  to  Problem  2:  given  a  pseudoflow  /,  find  the  e  such  that  /  is  e-tight.  We  need  a 
definition.  For  a  directed  graph  G  with  arc  cost  function  c,  the  minimum  cycle  mean  of  G,  denoted 
by  /x(G,c),  is  the  minimum,  over  all  cycles  T  in  G,  of  the  mean  cost  of  T,  defined  to  be  the  total 
arc  cost  c(r)  of  T  divided  by  the  number  of  arcs  it  contains.  The  connection  between  minimum 
cycle  means  and  tight  error  parameters  is  given  by  the  following  theorem,  which  was  discovered  by 
Engel  and  Schneider  [22]  and  later  by  Goldberg  and  Tarjan  [47]: 

Theorem  4.2.3  [22]  Suppose  a  pseudoflow  /  is  not  optimal.  Then  /  is  e-tight  for  e  =  -/i(G/,c). 

Proof:  Assume  /  is  not  optimal.  Consider  any  cycle  T  in  G/.  Let  the  length  of  T  be  l.  For  any  e, 
let  c(')  be  the  cost  function  defined  above:  c^\v,w)  =  c(v,tv)  +  e  for  (v,  w)  £  Ef.  Let  e  be  such 
that  /  is  e-tight,  and  let  /x  =  /x(G/,c).  By  Theorem  4.2.2,  0  <  c^)(T)  =  c(T)  +  /e,  t.e.,  c(T)/l  >  -e. 
Since  this  is  true  for  any  cycle  T,  /x  >  -e,  i.e.,  e  >  -p.  Conversely,  for  any  cycle  T,  c(T)//  >  p, 
which  implies  td~^(r)  >  0.  By  Theorem  4.2.2,  this  implies  -/x  >  e.  | 

Karp  [66]  observed  that  the  minimum  mean  cycle  can  be  computed  in  0{nm)  time  by  extracting 
information  from  a  single  run  of  the  Bellman-Ford  shortest  path  algorithm.  This  gives  the  fastest 
known  strongly  polynomial  algorithm  for  computing  the  minimum  cycle  mean.  The  fastest  scaling 
algorithm  is  the  0(>/ji”ilog(nC))-time  algorithm  of  Orlin  and  Ahuja  [84].  Since  we  are  interested 
here  in  strongly  polynomial  algorithms,  we  shall  use  Karp’s  bound  of  O(nm)  as  an  estimate  of  the 
time  to  compute  a  minimum  cycle  mean. 

The  following  observation  is  helpful  in  the  analysis  to  follow'.  Suppose  /  is  an  e-tight  pseudoflow 
and  c  >  0.  Let  p  be  a  price  function  such  that  /  is  c-optimal  with  respect  to  p,  and  let  T  be  a  cycle 
in  Gj  with  mean  cost  —  e.  Since  —  e  is  a  lower  bound  on  the  reduced  cost  of  an  arc  in  G/,  every  arc 
of  T  must  have  reduced  cost  exactly  —  e. 


4.3  Fixed  Arcs 

The  main  observation  that  leads  to  strongly  polynomial  cost-scaling  algorithms  for  tiie  minimum- 
cost  circulation  problem  is  the  following  result  of  Tardos  [95]:  if  the  absolute  value  of  the  reduced 
cost  of  an  arc  is  significantly  greater  then  the  current  error  parameter  e,  then  the  value  of  any 
optimal  circulation  on  this  arc  is  the  same  as  the  value  of  the  current  circulation.  The  following 
theorem  is  a  slight  generalization  of  this  result  (to  get  the  original  result,  take  e'  -  0).  This  theorem 
can  be  used  in  two  ways.  The  first  is  to  drop  the  capacity  constraint  for  an  arc  of  large  reduced 
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cost.  This  approach  is  used  in  [95].  The  second,  discussed  below,  is  to  consider  the  arcs  that  have 
the  same  flow  value  in  every  e-optimal  circulation  for  the  current  value  of  the  error  parameter  e 
and  to  notice  that  the  flow  through  these  arcs  will  not  change.  This  approach  is  used  in  [46,  43]. 

Theorem  4.3.1  [95]  Let  e  >  0,  d  >  0  be  constants.  Suppose  that  a  circulation  /  is  c-optimal  with 
respect  to  a  price  function  p,  and  that  there  is  an  arc  (v,w)  £  E  such  that  )ep(u,u;)|  >  n(e-fe').  Then, 
for  any  e'-optimal  circulation  /',  we  have  f(v,w)  =  f'(v,w). 

Proof:  By  antisymmetry,  it  is  enough  to  prove  the  theorem  for  the  case  cp(v,w)  >  n(e  +  d). 
Let  f  be  a  circulation  such  that  f(v,w)  ^  f(v,w).  Since  cp(v,w)  >  c,  the  flow  /  through  the 
arc  (v,w)  must  be  as  small  as  the  capacity  constraints  allow,  namely  —u(w,v),  and  therefore 
f'(v,w)  ^  f(v,w )  implies  f(v,w)  >  f(v,w).  We  show  that  /'  is  not  e'-optimal,  and  the  theorem 
follows. 

Consider  G>  =  {(x,j/)  G  E\f\x,y)  >  f{x,y)}.  Note  that  G>  is  a  subgraph  of  G/,  and  (t\  w) 
is  ai.  arc  of  G>.  Since  /  and  /'  are  circulations,  G>  must  contain  a  simple  cycle  T  that  passes 
through  (v,w).  Let  l  be  the  length  of  T.  Since  all  arcs  of  T  are  residual  arcs,  the  cost  of  T  is  at 
least 

cp(v,w)  -  (l  -  l)e  >  n(e  -|-  e')  -  (n  -  l)c  >  nc'. 

Now  consider  the  cycle  F  obtained  by  reversing  the  arcs  on  T.  Note  that  F  is  a  cycle  in 
G<  =  {(r,y)  G  E\f(x,y)  <  f(x,y)}  and  is  therefore  a  cycle  in  G/'.  By  antisymmetry,  the  cost  of 
T  is  less  than  —nP  and  thus  the  mean  cost  of  T  is  less  than  —  But  Theorem  4.2.3  implies  that 
/'  is  not  e'-optimal.  | 

To  state  an  important  corollary  of  Theorem  4.3.1,  we  need  the  following  definition.  We  say  that 
an  arc  (r,  w)  €  E  is  (-fixed  if  the  flow  through  this  arc  is  the  same  for  all  e-optimal  circulations. 

Corollary  4.3.2  [46]  Let  c  >  0,  suppose  /  is  e-optimal  with  respect  to  a  price  function  p,  and  suppose 
that  (v,w)  is  an  arc  such  that  |cp(u,u;)|  >  2 ne.  Then  ( v,w )  is  e-fixed. 

Define  F(  to  be  the  set  of  e-fixed  arcs.  Since  the  generalized  cost-scaling  method  decreases  e, 
an  arc  that  becomes  e-fixed  stays  e-fixed.  We  show  that  when  e  decreases  by  a  factor  of  2n,  a  new 
arc  becomes  e-fixed. 

Lemma  4.3.3  Assume  e'  <  Suppose  that  there  exists  an  e-tight  circulation  /.  Then  Fe<  properly 
contains  Ft. 

Proof:  Since  every  e'-optimal  circulation  is  e-optimal,  we  have  Ft  C  Ft>.  To  show  that  the  contain¬ 
ment  is  proper,  we  have  to  show  that  there  is  an  e'-fixed  arc  that  is  not  e-fixed. 
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Since  the  circulation  /  is  e-tight,  there  exists  a  price  function  p  such  that  /  is  c-optimal  with 
respect  to  p,  and  there  exists  a  simple  cycle  r  in  G/  every  arc  of  which  has  reduced  cost  — c.  (See 
Section  4.2.)  Since  increasing  /  along  V  preserves  £-optimality,  the  arcs  of  T  tire  not  c-fixed. 

We  show  that  at  least  one  arc  of  T  is  e'-fixed.  Let  /'  be  a  circulation  that  is  c'-optimal  with 
respect  to  some  price  function  p'.  Since  the  mean  cost  of  T  is  — e,  there  is  an  arc  {v,w)  of  T  with 
cp'(v,w )  <  —  e  <  —2n<J.  By  Corollary  4.3.2,  the  arc  (v,w)  is  e'-fixed.  | 

In  the  next  section  we  show  how  to  use  this  lemma  to  get  a  strongly  polynomial  bound  for  a 
variation  of  the  generalized  cost-scaling  method. 


4.4  The  Strongly  Polynomial  Framework 

The  minimum-cost  circulation  framework  of  Section  3.3  has  the  disadvantage  the*  the  number  of 
iterations  of  refine  depends  on  the  magnitudes  of  the  costs.  If  the  costs  are  huge  integer:,  *he 
method  need  not  run  in  time  polynomial  in  n  and  m\  if  the  costs  are  irrational,  the  method  need 
not  even  termina*e.  In  this  section  we  show  that  a  natural  modification  of  the  generalized  cost¬ 
scaling  approach  produces  strongly  polync...iol  algorithms.  The  running  time  bounds  we  derived  for 
algorithms  based  on  the  approach  of  Section  3.3  remain  valid  for  the  modified  approach  presented 
in  ili is  section.  The  main  idea  of  this  modification  is  to  improve  c  periodically  by  finding  a  price 
function  that  fits  the  current  circulation  better  than  the  current  price  function.  This  idea  can  also 
be  used  to  improve  the  practical  performance  of  the  method. 

The  changes  needed  to  make  the  generalized  cost-scaling  approach  strongly  polynomial,  sug¬ 
gested  by  Lemma  4.3.3,  are  to  add  an  extra  computation  to  the  main  loop  of  the  algorithm  and 
to  change  the  termination  condition.  Before  calling  refine  to  reduce  the  error  parameter  c,  the 
new  method  computes  the  value  A  and  a  price  function  p\  such  that  the  current  circulation  /  is 
A- tight  with  respect  to  p\.  The  strongly  polynomial  method  is  described  on  Figure  4.1.  The  value 
of  A  and  the  price  function  p\  in  line  (*)  are  computed  as  described  in  Section  4.2.  The  algorithm 
terminates  when  the  circulation  /  is  optimal,  i.e.,  A  =  0. 

The  time  to  perform  line  (*)  is  O(nm).  (See  Section  4.2.)  Since  all  the  implementations  of 
refine  that  we  have  considered  have  a  time  bound  greater  than  O(nm),  the  time  per  iteration  in 
the  new  version  of  the  algorithm  exceeds  the  time  per  iteration  in  the  original  version  by  less  than 
a  constant  factor.  Since  each  iteration  at  least  halves  e,  the  bound  of  0(log(nC))  on  the  number  of 
iterations  derived  in  Chapter  3  remains  valid,  assuming  that  the  costs  are  integral.  For  arbitrary 
real-valued  costs,  we  shall  derive  a  bound  of  0(m  log  n)  on  the  number  of  iterations. 

Theorem  4.4.1  {45]  The  total  number  of  iterations  of  the  while  loop  in  procedure  min-cost  is 
O(mlogn). 

Proof :  Consider  a  time  during  the  execution  of  the  algorithm.  During  the  next  0(log  n)  iterations. 
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procedure  mtn-cosi(V,  E,  u,  c); 

[initialization] 

( <—  C; 

Vu,  p(t-)  —  0; 

V(t>,tu)  G  E,  f(v,  u’)  «—  0; 

[loop] 

while  (  >  0  do  begin 

(*)  find  A  and  p\  such  that  /  is  A-tight  with  respect  to  px; 

if  A  >  0  then  (e,f,p)  <—  refine(\,  f,p\) 

else  return(/); 

end. 


Figure  4.1:  The  strongly  polynomial  algorithm. 

either  the  algorithm  terminates,  or  the  error  parameter  is  reduced  by  a  factor  of  2 n.  In  the  latter 
case,  Lemma  4.3.3  implies  that  an  arc  becomes  fixed.  If  all  arcs  become  fixed,  the  algorithm 
terminates  in  one  iteration  of  the  loop.  Therefore  the  total  number  of  iterations  is  O(mlogrc).  | 

The  best  strongly  polynomial  implementation  of  the  generalized  cost-scaling  method  [46],  based 
on  the  dynamic  tree  implementation  of  refine ,  runs  in  0(nm2  \og(n2/m)  log  n)  time. 


4.5  Cycle- Canceling  Algorithms 

The  ideas  discussed  in  Sections  4.2  -  4.4  are  quite  powerful.  In  this  section  we  use  these  ideas  to 
show  that  a  simple  cycle-canceling  algorithm  of  Klein  [70]  becomes  strongly  polynomia  r  a  careful 
choice  is  made  among  possible  cycles  to  cancel.  Klein’s  algorithm  consists  of  repeatedly  finding 
a  residual  cycle  of  negative  cost  and  sending  as  much  flow  as  possible  around  the  cycle.  This 
algorithm  can  run  for  an  exponential  number  of  iterations  if  the  capacities  and  costs  are  integers, 
and  it  need  not  terminate  if  the  capacities  are  irrational  [26].  Goldberg  and  Tarjan  [45]  showed 
that  if  a  cycle  with  the  minimum  mean  cost  is  canceled  at  each  iteration,  the  algorithm  becomes 
strongly  polynomial.  We  call  the  resulting  algorithm  the  minimum-mean  cycle-canceling  algorithm. 

The  minimum-mean  cycle-canceling  algorithm  is  closely  related  to  the  shortest  augmenting  path 
maximum  flow  algorithm  of  Edmonds  and  Karp  [21].  The  relationship  is  as  follows.  If  a  maximum 
flow  problem  is  formulated  as  a  minimum-cost  circulation  problem  in  a  standard  way,  then  Klein's 
cycle-canceling  algorithm  corresponds  exactly  to  the  Ford-Fulkerson  maximum  flow  algorithm,  and 
the  minimum-mean  cycle-canceling  algorithm  corresponds  exactly  to  the  Edmonds-Karp  algorithm. 
The  minimum-mean  cycle-canceling  algorithm  can  also  be  interpreted  as  a  steepest  descent  method 
using  the  L\  metric. 

For  a  circulation  /  we  define  c(f)  to  be  zero  if  /  is  optimal  and  to  be  the  unique  number  t'  >  0 
such  that  /  is  c'-  tight  otherwise.  We  use  e{f)  as  a  measure  of  the  quality  of  /.  Let  /  be  an  arbitrary 
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circulation,  let  e  =  ((f),  and  let  p  be  a  price  function  with  respect  to  which  /  is  c-optimal.  Holding 
e  and  p  fixed,  we  study  the  effect  on  e(f)  of  a  minimum-mean  cycle  cancellation  that  modifies  /. 
Since  all  arcs  on  a  minimum-mean  cycle  have  negative  reduced  cost  with  respect  to  p.  cancellation 
of  such  a  cycle  does  not  introduce  a  new  residual  arc  with  negative  reduced  cost,  and  hence  ((  f) 
does  not  increase. 

Lemma  4.5.1  A  sequence  of  m  minimum-mean  cycle  cancellations  reduces  ((f)  to  at  most  (1-1  /n)e , 
i.e.,  to  at  most  1  -  1/n  times  its  original  value. 

Proof:  Let  p  a  price  function  such  that  /  is  e-tight  with  respect  to  p.  Holding  e  and  p  fixed,  we 
study  the  effect  on  the  admissible  graph  G a  (with  respect  to  the  circulation  /  and  price  function  p) 
of  a  sequence  of  m  minimum-mean  cycle  cancellations  that  modify  /.  Initially  every  arc  (v,  tv)  £  £4 
satisfies  cp(v,  w)  >  — c.  Canceling  a  cycle  all  of  whose  arcs  are  in  Ea  adds  only  arcs  of  positive 
reduced  cost  to  Ej  and  deletes  at  least  one  arc  from  E a-  We  consider  two  cases. 

Case  1:  None  of  the  cycles  canceled  contains  an  arc  of  nonnegative  reduced  cost.  Then  each 
cancellation  reduces  the  size  of  Ea,  and  after  m  cancellations  £4  is  empty,  which  implies  that  /  is 
optimal,  i.e.,  ((f)  =  0.  Thus  the  lemma  is  true  in  this  case. 

Case  2:  Some  cycle  canceled  contains  an  arc  of  nonnegative  reduced  cost.  Let  T  be  the  first 
such  cycle  canceled.  Every  arc  of  T  has  a  reduced  cost  of  at  least  -c,  one  arc  of  V  has  a  nonnegative 
reduced  cost,  and  the  number  of  arcs  in  T  is  at  most  n.  Therefore  the  mean  cost  of  T  is  at  least 
—  (1  —  l/n)e.  Thus,  just  before  the  cancellation  of  T,  ((f)  <  (1  —  l/n)e  bv  Theorem  4.2.3.  Since 
((f)  never  increases,  the  lemma  is  true  in  this  case  as  well.  | 

Lemma  4.5.1  is  enough  to  derive  a  polynomial  bound  on  the  number  of  iterations,  assuming 
that  all  arc  costs  are  integers. 

Theorem  4.5.2  [45]  If  all  arc  costs  are  integers,  then  the  minimum-mean  cycle-canceling  algorithm 
terminates  after  0(nm  log(nC))  iterations. 

Proof:  The  lemma  follows  from  Lemmas  3.2.1  and  4.5.1  and  the  observation  that  the  initial  circu¬ 
lation  is  C-optimal.  | 

To  obtain  a  strongly  polynomial  bound,  we  use  the  ideas  of  Section  4.3.  The  proof  of  the  next 
theorem  uses  the  following  inequality: 

(1  -  ^)n(lnn'M)  <  ±  for  n  >  2. 

Theorem  4.5.3  [45]  For  arbitrary  real-valued  arc  costs,  the  minimum-mean  cycle-canceling  algorithm 
terminates  after  0(nm2  log  n)  cycle  cancellations. 
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Proof:  Let  k  =  m(n  fin  n  +  1]),  Divide  the  iterations  into  groups  of  k  consecutive  iterations.  We 
claim  that  each  group  of  iterations  fixes  the  flow  on  a  distinct  arc  (v,w),  i.e .,  iterations  after  those 
in  the  group  do  not  change  f(v,  w).  The  theorem  is  immediate  from  the  claim. 

To  prove  the  claim,  consider  any  group  of  iterations.  Let  /  be  the  flow  before  the  first  iteration 
of  the  group,  f  the  flow  after  the  last  iteration  of  the  group,  t  —  ((f),  ('  =  ((f),  and  let  p'  be  a 
price  function  for  which  f  satisfies  the  e'-optimality  constraints.  Let  T  be  the  cycle  canceled  in  the 
first  iteration  of  the  group.  By  Lemma  4.5.1,  the  choice  of  k  implies  that  ('  <  c(l  -  i)nrini>+il  < 
4^.  Since  the  mean  cost  of  T  is  -e,  some  arc  on  F,  say  (u,tn),  must  have  cv>(v,w)  <  —  e  <  -2 nP . 
By  Corollary  4.3.2,  the  flow  on  (v,w)  will  not  be  changed  by  iterations  after  those  in  the  group. 
But  f(v,  u>)  is  changed  by  the  first  iteration  in  the  group,  which  cancels  T.  Thus  each  group  fixes 
the  flow  on  a  distinct  arc.  | 


Theorem  4.5.4  [45]  The  minimum-mean  cycle-canceling  algorithm  runs  in  0(n2m3Iogn)  time  on 
networks  with  arbitrary  real-valued  arc  costs,  and  in  0(n2m2  min{log(rjC).  m  log  n})  time  on  networks 
with  integer  arc  costs. 

Proof :  Immediate  from  Theorems  4.5.2  and  4.5.3.  | 

Although  the  minimum-mean  cycle-canceling  algorithm  seems  to  be  of  mostly  theoretical  in¬ 
terest,  it  has  a  variant  that  is  quite  efficient.  This  variant  maintains  a  price  function  and.  instead 
of  canceling  the  minimum-mean-cost  residual  cycle,  cancels  residual  cycles  composed  entirely  of 
negative  reduced  cost  arcs;  if  no  such  cycle  exists,  the  algorithm  updates  the  price  function  to 
improve  the  error  parameter  c.  An  implementation  of  the  algorithm  using  the  dynamic  tree  data 
structure  runs  in  0(nm  log  n  log(nC))  time.  See  [45]  for  details. 

The  techniques  used  to  analyze  the  minimum-mean  cycle-canceling  algorithm  also  provide  an 
approach  to  making  the  primal  network  simplex  algorithm  for  the  minimum-cost  circulation  prob¬ 
lem  more  efficient.  Tarjan  [99]  has  discovered  a  pivot  rule  that  gives  a  bound  of  0(n*  Iosr,+<9(|)) 
on  the  number  of  pivots  and  the  total  running  time.  This  is  the  first  known  subexponential  time 
bound.  For  the  extended  primal  network  simplex  method  in  which  cost-increasing  as  well  as  cost- 
decreasing  pivots  are  allowed,  he  has  obtained  a  strongly  polynomial  time  bound.  In  contrast, 
the  dual  network  simplex  method  is  known  to  have  strongly  polynomial  variants,  as  mentioned  in 
Chapter  5. 

Barahona  and  Tardos  [7]  have  exhibited  another  cycle-canceling  algorithm  that  runs  in  polyno¬ 
mial  time.  Their  algorithm  is  based  on  an  algorithm  of  Weintraub  [103],  W'hich  works  as  follows. 
Consider  the  improvement  in  the  cost  of  a  circulation  obtained  by  canceling  a  negative  cycle.  Since 
a  symmetric  difference  of  two  circulations  can  be  decomposed  into  at  most  m  cycles,  canceling  the 
cycle  that  gives  the  best  improvement  reduces  the  difference  between  the  current  and  the  optimal 
values  of  the  cost  by  a  factor  of  (1  -  1/m).  If  the  input  data  is  integral,  only  a  polynomial  number 
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of  such  improvements  can  be  made  until  an  optimal  solution  is  obtained.  Finding  the  cycle  that 
gives  the  best  improvement  is  NP-hard,  however.  Weintraub  shows  how  to  find  a  collection  of 
cycles  whose  cancellation  reduces  the  cost  by  at  least  as  much  as  the  best  improvement  achievable 
by  canceling  a  single  cycle.  His  method  requires  a  superpolynomial  number  of  applications  of  an 
algorithm  for  the  assignment  problem.  Each  such  application  yields  a  minimum-cost  collection  of 
vertex-disjoint  cycles.  Barahona  and  Tardos  [7]  have  shown  that  the  algorithm  can  be  modified  so 
that  the  required  collection  of  cycles  is  found  in  at  most  m  assignmen*  computations.  The  resulting 
minimum-''ost  circulation  algorithm  runs  in  polynomial  time. 
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Chapter  5 


Capacity- Scaling  Algorithms 


5.1  Introduction 


In  this  chapter,  we  survey  minimum-cost  circulation  algorithms  based  on  capacity  scaling.  We 
concentrate  on  two  algorithms:  the  first  polynomial-time  algorithm,  that  of  Edmonds  and  Karp  [21], 
who  introduced  the  idea  of  scaling,  and  the  algorithm  of  Orlin  [82].  The  latter  is  an  extension  of 
the  Edmonds-Karp  algorithm  and  is  the  fastest  known  strongly  polynomial  algorithm. 

In  this  chapter  we  shall  consider  the  (uncapacitated)  transshipment  problem,  which  is  equivalent 
to  the  minimum-cost  circulation  problem.  This  simplifies  the  presentation  considerably.  We  begin 
by  describing  a  generic  augmenting  path  algorithm,  which  we  call  the  minimum-cost  augmentation 
algorithm,  that  is  the  basis  of  most  capacity-scaling  algorithms.  The  algorithm  is  due  to  Jewel  [61], 
Busacker  and  Gowen  [14],  and  Iri  [58].  The  use  of  dual  variables  as  described  below  was  proposed 
independently  by  Edmonds  and  Karp  [21]  and  Tomizawa  [100].  The  idea  of  the  algorithm  is  to 
maintain  a  pseudoflow  /  that  satisfies  the  complementary  slackness  constraints  while  repeatedly 
augmenting  flow  to  gradually  get  rid  of  all  excesses.  Two  observations  justify  this  method:  (1) 
augmenting  flow  along  a  minimum-cost  path  preserves  the  invariant  that  the  current  pseudoflow 
has  minimum  cost  among  those  pseudoflows  with  the  same  excess  function;  (2)  a  shortest  path 
computation  suffices  both  to  find  a  path  along  which  to  move  flow  and  to  find  price  changes  to 
preserve  complementary  slackness. 

The  algorithm  maintains  a  pseudoflow  /  and  a  price  function  p  such  that  cp(v,w)  >  0  for 
every  residual  arc  (v,w).  The  algorithm  consists  of  repeating  the  augmentation  step ,  described  in 
Figure  5.1,  until  every  excess  is  zero.  Then  the  current  pseudofiow  is  optimal. 

The  correctness  of  this  algorithm  follows  from  the  observation  that  the  price  transformation 
in  the  augment  step  preserves  complementary  slackness  and  makes  the  new  reduced  cost  of  any 
minimum-cost  path  from  s  to  t  equal  to  zero.  One  iteration  of  the  augment  step  takes  time 
proportional  to  that  required  by  a  single-source  shortest  path  computation  on  a  graph  with  arcs 
of  non-negative  cost.  The  fastest  known  strongly  polynomial  algorithm  for  this  computation  is 
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augments,  t). 

Applicability:  e(s)  >  0  and  e(t)  <  0. 

Action:  For  every  vertex  v,  compute  rr(v),  the  minimum  reduced  cost  of  a  residual 
path  from  s  to  v.  For  every  vertex  v,  replace  p(v)  by  p(v)  +  ?r(t>).  Move  a  positive 
amount  of  flow  from  s  to  t  along  a  path  of  minimum  reduced  cost. 


Figure  5.1:  The  augmentation  step. 

Fredman  and  Tarjan’s  implementation  [29]  of  Dijkstra’s  algorithm,  which  runs  in  0(nlogn  +  m) 
time.  The  problem  can  also  be  solved  in  O(mloglogC)  time  [85]  or  0(ny/\ogC  +  m)  time  [2], 
where  C  is  the  maximum  arc  cost.  Since  we  are  mainly  interested  here  in  strongly  polynomial 
algorithms,  we  shall  use  0(nlogn  +  m)  as  our  estimate  of  the  time  for  each  augment  step. 

Remark:  The  only  reason  to  maintain  prices  in  this  algorithm  is  to  simplify  the  minimum-cost  path 
computations  by  guaranteeing  that  each  such  computation  is  done  on  a  graph  all  of  whose  arc  costs 
are  non-negative.  If  prices  are  not  maintained  in  this  way,  each  shortest  path  computation  takes 
O(nm)  time  using  the  Bellman-Ford  algorithm  (see  e.g.  [96]). 

5.2  The  Edmonds-Karp  Algorithm 

The  overall  running  time  of  the  augmentation  algorithm  depends  on  the  number  of  augmentation 
steps,  which  cannot  be  polynomially  bounded  without  specifying  their  order  more  precisely.  To  im¬ 
pose  an  efficient  order,  Edmonds  and  Karp  introduced  the  idea  of  capacity  scaling,  which  Orlin  [80] 
in  his  description  of  Edmonds-Karp  algorithm  reinterpreted  as  excess  scaling.  We  shall  present  a 
modification  of  Orlin’s  version  of  the  Edmonds-Karp  algorithm. 

The  algorithm  maintains  a  scaling  parameter  A  such  that  the  flow  on  every  arc  is  an  integer 
multiple  of  A.  For  a  given  pseudoflow  /  and  value  of  the  scaling  parameter  A,  we  denote  the  sets 
of  vertices  with  large  excesses  and  large  deficits  as  follows: 

57(A)  =  {v  6  V  :  ef(v)  >  A}; 

77(A)  =  {i>  G  E  :  e]{v)  <  -A}.  (5.1) 

The  algorithm  consists  of  a  number  of  scaling  phases.  During  a  A-scaling  phase  the  residual 
capacity  of  every  arc  is  an  integer  multiple  of  A.  The  algorithm  chooses  vertices  s  G  5/(  A/2),  1  G 
Tj(  A/2),  performs  augments,  t),  which  sends  A  units  of  flow  along  a  minimum-cost  path  from  s 
to  t,  and  repeats  such  augmentations  until  either  S/( A/2)  or  7/(A/2)  is  empty.  Note  that  this 

can  create  a  deficit  in  place  of  an  excess,  and  vice  versa.  Vertices  with  the  new  excesses  and 

deficits,  however,  are  not  in  S/(A/2)  \JTj(A/2).  Then  the  algorithm  halves  A  and  begins  a  new 

scaling  phase.  The  initial  value  of  A  is  2^ogDl.  The  algorithm  terminates  after  the  1-scaling  phase, 
assuming  all  supplies  are  integral.  Figure  5.2  provides  a  detailed  description  of  a  phase  of  the 
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while  5, (A/2)  ^  0  and  7}(A/2)  ^  0  do  begin 

Choose  s  £  S/(A/2)  and  t  £  Tj( A/2).  Perform  augment  ( s,t ),  sending  A  units  of 
flow  along  the  augmenting  path  selected. 

end; 


I 


Figure  5.2:  A  phase  of  the  Edmonds-Karp  Algorithm. 


algorithm. 

The  following  lemma  follows  from  the  description  of  the  algorithm. 

Lemma  5.2.1  During  a  A-scaling  phase,  the  residual  capacity  of  every  arc  is  an  integer  multiple  of 
A. 


Lemma  5.2.2  Let  /'  be  the  pseudoflow  at  the  end  of  a  2A-scaling  phase.  If  S/<(A)  =  0,  then  there 
are  at  most  |S/<(A/2)|  flow  augmentations  during  the  A-scaling  phase.  An  analogous  statement  is  true 
if  Tf(  A)  =  0  at  the  end  of  a  2A-scaling  phase. 

Proof:  By  assumption,  Sj>{ A)  =  0.  Therefore,  for  every  vertex  v  £  V,  ej( v)  <  A  during  the 
A-scaling  phase.  This  implies  that  pushing  A  units  of  flow  along  a  path  from  s  £  5/(A/2)  to 
t  £  T/(A/2)  removes  s  from  Sj( A/2).  Note  that  since  e/(v)  <  -A/2  for  v  £  T( A/2),  vertices  with 
new  excesses  are  not  in  S/( A/2).  | 

There  are  [logD]  phases,  each  of  which  consists  of  up  to  n  single-source  shortest  path  compu¬ 
tations  on  networks  with  non-negative  arc  costs.  Thus  we  have  the  following  results  for  networks 
with  integral  supplies: 


Theorem  5.2.3  [21]  The  Edmonds-Karp  algorithm  solves  the  transshipment  problem  in 
0((n  log  D)  (n  log  n  +  m))  time,  and  the  minimum-cost  circulation  problem  in  0((m  log  U)(n  log  n  + 
m))  time. 


The  excess-scaling  idea  can  be  combined  with  the  idea  of  augmenting  along  approximately 
minimum-cost  paths  rather  than  exactly  minimum-cost  paths.  Approximately  minimum-cost  paths 
can  be  defined  using  the  c-optimality  notion  discussed  in  Chapter  3.  An  appropriate  combination  of 
these  ideas  yields  the  double  sealing  algorithm  of  Ahuja,  Goldberg,  Orlin,  and  Tarjan  [1].  An  imple¬ 
mentation  of  this  algorithm  based  on  dynamic  trees  has  a  time  bound  of  O(nm(log log  V)  log(nC)) 
for  the  minimum-cost  circulation  problem. 
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5.3  Strongly  Polynomial  Algorithms 

After  Tardos  [95]  discovered  the  first  strongly  polynomial  algorithm,  Fujishige  [30]  and  indepen¬ 
dently  Orlin  [80]  developed  more  efficient  strongly  polynomial  algorithms.  Orlin’s  algorithm  [82], 
that  is  the  main  focus  of  this  section,  is  an  extension  of  these  algorithms.  These  algorithms  are 
based  on  the  method  of  Section  5.2  and  the  following  dual  version  of  Theorem  4.3.1. 

Theorem  5.3.1  [30,  80]  Let  /  be  a  pseudoflow,  and  let  p  be  a  price  function  such  that  the  comple¬ 
mentary  slackness  conditions  are  satisfied  and  f(v,w)  >  Hv:e/(„)>o  le/(u)l-  Then  every  optimal  price 
function  p*  satisfies  cp»(v,w)  =  0. 

Proof:  Let  p*  be  an  optimal  price  function,  and  let  /*  be  a  corresponding  optimal  pseudofiow. 
Assume  by  way  of  contradiction  that  cp*( w,  v)  ^  0.  Define  the  pseudoflow  f  by  f'(x,  y)  =  f(x ,  y )  — 
f’(x,y)  (i.e.,  we  can  obtain  an  optimal  pseudoflow  by  augmenting  /  by  /').  The  Decomposition 
Theorem  (Theorem  1.7.2)  i^nlies  that  /'  r* n  he  d^omposed  into  flows  along  a  collection  of  cycles 
and  a  collection  of  paths  from  vertices  with  excess  to  vertices  with  deficits  (both  excesses  and 
deficits  are  with  respect  to  /).  Furthermore,  the  Decomposition  Theorem  also  implies  that  for  any 
arc  (x,y)  that  is  not  on  one  of  the  cycles,  |/'(x,y)|  <  E»:e;(„)>o  ef(v)- 

The  arcs  on  the  cycles  of  the  decomposition  are  in  Ej,  and  the  arcs  opposite  to  the  ones  on 
the  cycles  are  in  Ef.  This  implies  that  the  cycles  have  zero  cost,  and  that  every  arc  on  the  cycles 
must  have  zero  reduced  cost  with  respect  to  both  p  and  p*.  We  have  assumed  that  cp.(u,  w)  0; 
therefore  (w,v)  cannot  be  on  such  a  cycle.  Thus  f(w,v)  <  ]£v:e/(v)>o  ef(v)  and  thus  /*(t\u’)  >  0. 
But  f‘(v,w)  >  0  and  cp»(v,w)  /  0  contradict  the  complementary  slackness  constraint.  | 

The  idea  behind  the  strongly  polynomial  algorithms  based  on  capacity-scaling  is  to  contract  the 
arcs  satisfying  the  condition  o:  Theorem  5.3.1.  Let  /  and  p  be  a  pseudoflow  and  a  price  function, 
respectively,  and  let  {v,w)  be  an  arc  satisfying  the  condition  of  Theorem  5.3-1.  The  reduced  cost 
function  cp  defines  an  equivalent  problem  with  cp(v,w)  =  0.  By  the  above  theorem  the  optimal 
prices  of  the  vertices  v  and  w  are  the  same.  Define  a  transshipment  problem  on  the  network  formed 
by  contracting  the  arc  (u,ic),  with  cost  function  cp,  capacity  function  u,  and  demand  function  d 
such  that  the  demand  of  the  new  vertex  is  d(v)  +  d(w). 

Theorem  5.3.2  For  any  optimal  price  function  p*  for  the  new  problem,  the  price  function  p'(r)  = 
p(r )  +  p*(r )  for  r  G  V  (with  p*(i>)  =  p*(w)  defined  to  be  the  price  of  the  new  vertex)  is  optimal  for 
the  original  problem. 

Proof:  The  optimal  price  function  is  the  optimal  solution  of  a  linear  program,  the  linear  pro¬ 
gramming  dual  of  the  minimum-cost  flow  problem.  The  price  function  p'  =  p  -f  p*  is  the  optimal 
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Step  1  Run  the  Edmonds-Karp  algorithm  for  the  first  [21ogn]  scaling  phases  with  cost 
function  Cp,  starting  with  A  =  max^y  |d(t>)|.  Let  /'  be  the  pseudoflow  and  let  p' 
be  the  price  function  found  by  the  algorithm. 

Step  2  Contract  every  arc  ( v,w )  with  \f(v,w)\  >  nA  and  update  the  price  function  by 
setting  p(v)  <—  p(v)  +  p'(v )  for  all  v  6  V  (where  all  vertices  v  €  V  contracted  into 
the  same  vertex  t/  of  the  current  network  have  the  same  price  p'{v)  =  p'(v')). 


Figure  5.3:  The  inner  loop  of  the  simple  strongly  polynomial  algorithm. 

solution  of  same  linear  program  with  the  additional  constraint  c(v,  w)  +  p'(v)  -  p'{w )  =  0.  By 
Theorem  5.3.1,  the  two  linear  programs  have  the  same  set  of  optimal  solutions.  | 

We  shall  describe  a  simple  strongly  polynomial  algorithm  based  on  this  idea.  It  is  a  variation 
of  Fujishige’s  algorithm  [30],  though  our  description  is  simplified  by  using  ideas  from  [82].  One 
iteration  of  the  algorithm  consists  of  running  the  Edmonds-Karp  algorithm  for  the  first  21ogn 
scaling  phases  starting  with  A  =  max^gv'  | cf ( v ) |  (we  cannot  find  the  power  of  two  used  by  the 
Edmonds-Karp  algorithm  in  strongly  polynomial  time),  and  then  contracting  all  arcs  that  satisfy 
the  condition  of  Theorem  5.3.1.  The  next  iteration  considers  the  contracted  network.  We  shall 
prove  that  °ach  iteration  will  contract  at  least  one  arc.  The  algorithm  terminates  when  the  current 
pseudoflow  is  a  feasible  solution.  Theorem  5.3.2  guarantees  that  the  price  function  found  at  this 
point  is  optimal.  Given  an  optimal  price  function,  the  optimal  flow  can  be  found  by  a  single 
maximum  flow  computation.  (See  Figure  5.3  for  a  more  detailed  description.) 

The  following  theorem  implies  that  at  least  one  arc  is  contracted  at  each  iteration  of  the  inner 
loop. 

Theorem  5.3.3  [82]  Let  /  and  A  be  the  pseudoflow  and  the  scaling  parameter  at  the  end  of  an 
iteration  of  the  algorithm  of  Figure  5.3.  Then  le/(l’)l  <  nA.  Furthermore,  if  /  is  not  feasible, 

then  there  exists  an  arc  (r,  re)  with  \f(v,  ic)|  >  nA. 

Proof:  After  the  A-scaling  phase  either  S/( A/2)  or  Tj( A/2)  is  empty.  For  a  pseudoflow  /,  the 
excesses  sum  to  zero.  Therefore  |®/(n)|  <  nA.  Also,  every  excess  is  less  than  nA.  After 

[21ogn]  scaling  phases,  A  <  (maxwev  |d(w)|)/n2.  For  a  vertex  v  whose  demand  has  the  maximum 
absolute  value  this  implies  that 

I  5Z  f(v'w)\  ^  lrf(u)l  -nA>  n(n  -  1)A.  (5.2) 

wdV 


Consequently,  at  least  one  arc  incident  to  v  carries  at  least  nA  flow.  | 

One  iteration  takes  O(n(n\ogn  +  m)logn)  time.  Each  iteration  contracts  at  least  one  arc,  and 
therefore  there  are  at  most  n  iterations. 
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Theorem  5.3.4  The  algorithm  of  Figure  5.3  solves  the  transshipment  problem  in  0(n2  log  n(n  log  n  + 
m ))  time,  and  the  minimum-cost  circulation  problem  in  0(m2logn(nlogn  +  m))  time. 

This  simple  algorithm  wastes  a  lot  of  time  by  throwing  away  the  current  flow  after  each  con¬ 
traction.  If  there  are  several  vertices  with  demand  close  to  the  maximum  then  it  might  make  sense 
to  run  somewhat  more  than  21ogn  scaling  phases.  The  proof  of  Theorem  5.3.3  would  then  apply 
to  more  than  one  vertex,  and  more  than  one  arc  could  be  contracted  in  an  iteration.  This  idea 
does  not  help  much  if  there  are  not  enough  vertices  with  close-to-maximum  demand.  On  the  other 
hand,  by  Lemma  5.2.2  the  number  of  shortest  path  computations  during  a  phase  can  be  bounded 
in  terms  of  the  number  of  vertices  with  relatively  large  demand  (it  is  bounded  either  by  |S(A/2)| 
or  by  |T(A/2)|).  One  could  try  to  find  a  point  for  stopping  the  scaling  algorithm  that  balances  the 
number  of  shortest  path  computations  done  and  the  number  of  arcs  contracted. 

Galil  and  Tardos  [37]  developed  an  0(n2(nlog  +  m)  log  n)-time  minimum-cost  circulation  algo¬ 
rithm  based  on  this  idea.  In  fact,  the  balancing  technique  of  Galil  and  Tardos  [37],  together  with 
Orlin’s  [82]  proof  of  Theorem  5.3.3,  can  be  used  to  give  an  0(n(n  log  n  +  m)  log  n)- time  algorithm 
for  the  transshipment  problem.  Instead  of  pursuing  this  approach  further,  however,  we  shall  present 
a  variant  of  a  much  simpler  algorithm  due  to  Orlin  [82]  that  has  same  efficiency  and  does  not  rely 
on  delicate  balancing.  This  algorithm  contracts  arcs  during  the  scaling  phases,  without  starting  the 
scaling  from  scratch  after  each  contraction;  the  algorithm  finds  an  optimal  flow  directly,  without  us¬ 
ing  an  optimal  price  function.  The  algorithm  is  a  modification  of  the  Edmonds-Karp  algorithm,  but 
runs  in  O(nlogn)  iterations  instead  of  the  0(n  log  D)  iterations  of  the  Edmonds-Karp  algorithm. 

Orlin’s  algorithm  runs  the  Edmonds-Karp  scaling  algorithm  starting  with  A  =  max,,€v  |d(t’)l 
and  contracts  all  arcs  that  carry  at  least  4nA  flow.  The  scaling  algorithm  continues  as  long  as 
the  current  pseudoflow  is  non-zero  on  some  uncontracted  arc.  When  the  pseudoflow  is  zero  on 
every  arc,  the  scaling  is  restarted  by  setting  A  to  be  the  maximum  absolute  value  of  a  demand. 
Roughly  speaking,  a  vertex  v  is  the  start  of  at  most  21ogrc  shortest  path  computations  during  the 
algorithm.  For  a  vertex  v  of  the  original  graph,  this  follows  from  the  fact  that  v  £  Sj( A)  U  T;{ A) 
unless  |d(r)|  >  A,  and  therefore  an  arc  incident  to  v  will  be  contracted  at  most  21ogn  scaling 
phases  after  v  first  served  as  a  starting  vertex.  (See  the  proof  of  Theorem  5.3.3.)  This  may  not  be 
true  for  the  contracted  vertices,  however.  A  vertex  created  by  a  contraction  may  have  very  small 
demand  relative  to  its  current  excess. 

We  say  that  a  vertex  v  is  active  if  v  E  5/(A/2)  U  7/(A/2).  A  vertex  v  is  activated  by  a  A- 
scaling  phase  if  v  is  not  active  at  the  end  of  the  2A-scaling  phase,  and  it  becomes  active  during 
the  A-scaling  phase  (i.e.,  either  it  becomes  active  at  the  beginning  of  the  phase  or  the  vertex  is 
created  by  contraction  during  the  phase).  We  can  prove  the  following  lemma  using  a  proof  similar 
to  that  of  Lemma  5.2.2. 

Lemma  5.3.5  The  number  of  shortest  path  computations  during  a  phase  is  bounded  by  the  number 
of  vertices  activated  by  the  phase. 


5.3.  STRONGLY  POLYNOMIAL  ALGORITHMS 


61 


Step  1.  Let  A  =  max^y  |d(t>)|.  If  A  =  0  then  STOP. 

Step  2.  while  Sj( A/2)  ^  0  and  Tj( A/2)  0  do  begin 

Choose  s  €  Sj( A/2)  and  t  €  T/( A/2).  Perform  augments,  <),  sending  A  units  of 
flow  along  the  augmenting  path  selected.  Contract  every  arc  (r,u»)  with  f(v,w)  > 

4nA. 

end. 

Step  3.  If  /  is  zero  on  all  uncontracted  arcs  go  to  Step  1,  otherwise  let  A  =  A/2  and  go  to 
_ Step  2. _ 

Figure  5.4:  Orlin’s  Algorithm. 

Theorem  5.3.6  [82]  A  vertex  can  be  activated  at  most  [21ogn]  +  1  times  before  it  is  contracted. 

Proof:  A  vertex  can  be  activated  once  due  to  contraction.  When  a  vertex  v  is  activated  for  the 
second  time,  it  must  already  exist  at  the  end  of  the  previous  scaling  phase,  when  it  was  not  active. 
Let  A  be  the  scaling  parameter  in  the  phase  when  v  is  activated  for  the  second  time.  At  the 
beginning  of  this  phase,  A/2  <  |ey(u)|  <  A.  But  d(v)  -  ej(v)  is  an  integer  multiple  of  2A.  This 
implies  that  A/2  <  |ey(u)|  <  |d(u)|.  After  O(logn)  more  scaling  phases  the  scaling  parameter 
will  be  less  than  \d(v)\/4n2.  Using  an  argument  analogous  to  the  proof  of  Theorem  5.3.3  we  can 
conclude  that  some  arc  incident  to  v  is  contracted.  | 

The  previous  lemma  and  theorem  bound  the  number  of  shortest  path  computations  during  the 
algorithm.  All  other  work  done  is  linear  per  scaling  phase.  At  least  one  arc  is  contracted  in  each 
group  of  O(logn)  scaling  phases.  Therefore,  there  are  at  most  O(nlogn)  scaling  phases. 

Theorem  5.3.7  [82]  Orlin’s  algorithm  solves  the  transshipment  problem  in  0(n  logminjn,  D}(n  log  n  + 
m))  time,  and  the  minimum-cost  circulation  problem  in  0(m  log  min{n,  U}{n  log  n  +  m))  time. 

Remark:  In  contrast  to  the  simpler  strongly  polynomial  algorithm  discussed  earlier,  Orlin’s  al¬ 
gorithm  constructs  an  optimal  pseudoflow  directly,  without  first  constructing  an  optimal  price 
function.  To  see  this,  consider  the  amount  of  flow  moved  during  the  A-scaling  phase  for  some  value 
A.  Lemma  5.3.5  gives  a  2nA  bound.  Suppose  an  arc  is  contracted  during  the  A-scaling  phase. 
The  overall  amount  of  flow  that  is  moved  after  this  contraction  can  be  bounded  by  4nA.  Therefore 
the  pseudoflow  constructed  by  the  algorithm  is  feasible  in  the  uncontracted  network. 

Orlin  [80]  has  observed  that  capacity  scaling  ideas  can  be  used  to  guide  pivot  selection  in 
the  dual  network  simplex  algorithm  for  the  transshipment  problem.  More  recently  Orlin,  Plotkin 
and  Tardos  (personal  communication,  1989)  have  obtained  an  O(nm\og  min(n,  D))  bound  on  the 
number  of  pivots  based  on  such  a  pivot  selection  strategy. 


Chapter  6 

The  Generalized  Flow  Problem 

6.1  Introduction 


The  generalized  flow  problem  models  the  following  situation  in  financial  analysis.  An  investor  wants 
to  take  advantage  of  the  discrepancies  in  the  prices  of  securities  on  different  stock  exchanges  and  of 
currency  conversion  rates.  His  objective  is  to  maximize  his  profit  by  trading  on  different  exchanges 
and  by  converting  currencies.  The  generalized  flow  problem,  considered  in  this  chapter,  models  the 
above  situation,  assuming  that  a  bounded  amount  of  money  is  available  to  the  investor  and  that 
bounded  amounts  of  securities  can  be  traded  without  affecting  the  prices.  Vertices  of  the  network 
correspond  to  different  currencies  and  securities,  and  arcs  correspond  to  possible  transactions. 

The  generalized  flow  problem  was  first  considered  by  Jewell  [61].  This  problem  is  very  similar  to 
the  minimum-cost  circulation  problem,  and  several  of  the  early  minimum-cost  circulation  algorithms 
have  been  adapted  to  this  problem.  The  first  simple  combinatorial  algorithms  were  developed  by 
Jewell  [61]  and  Onaga  [79].  These  algorithms  are  not  even  pseudopolynomial,  however,  and  for 
real-valued  data  they  need  not  terminate.  Several  variations  suggested  in  the  early  70's  result  in 
algorithms  running  in  finite  (but  exponential)  time.  The  paper  by  Truemper  [101]  contains  an 
excellent  summary  of  these  results. 

The  generalized  flow  problem  is  a  special  case  of  linear  programming,  and  therefore  it  can  be 
solved  in  polynomial  time  using  any  general-purpose  linear  programming  algorithm,  such  as  the 
ellipsoid  method  [69]  and  the  interior-point  algorithms  [65,  88].  Kapoor  and  Vaidya  [64]  developed 
techniques  to  use  the  structure  of  the  matrices  that  arise  in  linear  programming  formulations 
of  network  flow  problems  to  speed  up  interior-point  algorithms.  The  first  two  polynomial-time 
combinatorial  algorithms  were  developed  by  Goldberg,  Plotkin,  and  Tardos  [42].  We  shall  review- 
one  of  their  algorithms  in  detail  and  sketch  the  other  one. 

The  current  fastest  algorithm  for  the  generali  d  flow  problem,  due  to  Vaidya,  is  based  on 
an  interior-point  method  for  linear  programming  and  runs  in  0(n2mh5  log(nB))  time  [102].  This 
algorithm  combines  the  ideas  from  the  paper  of  Kapoor  and  Vaidya  [64]  with  the  current  fastest 
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linear  programming  algorithm  (which  uses  fast  matrix  multiplication).  The  two  combinatorial 
algorithms  due  to  Goldberg,  Plotkin  and  Tardos  142]  run  in  0(mn2(m  +  nlogn)logrclogi?)  and 
0(m2n2  logo  log2  B)  time,  respectively.  More  recently  Maggs  (personal  communication,  1989) 
improved  the  latter  bound  through  better  balancing  to  0  ^n2m2/og2i?log^nslog^^^_(_log y— g ) . 

6.2  Simple  Combinatorial  Algorithms 

Jewell  [61]  and  Onaga  [78]  suggested  solving  the  generalized  flow  problem  by  using  adaptations  of 
augmenting  path  algorithms  for  the  maximum  flow  and  minimum-cost  circulation  problems.  The¬ 
orem  1.6.2  is  the  basis  of  Onaga’s  augmenting  path  algorithm  for  the  restricted  problem.  Starting 
from  the  zero  flow,  this  algorithm  iteratively  augments  the  flow  along  a  flow-generating  cycle  in  the 
residual  graph,  thereby  increasing  the  excess  at  the  source.  The  structure  of  the  restricted  problem 
makes  it  possible  to  find  the  most  efficient  flow-generating  cycle,  i.e.,  the  residual  flow-generating 
cycle  with  highest-gain.  The  highest-gain  cycle  consists  of  an  arc  (r, s)  and  a  highest-gain  simple 
path  from  s  to  r  for  some  vertex  r.  The  highest-gain  simple  path  to  r  is  the  shortest  path,  if  the 
length  of  each  arc  (i>,  w)  is  defined  as  l(v,w)  =  -log 7 (v,w).  Such  a  shortest  path  exists  since 
deleting  the  arcs  entering  the  source  from  the  residual  graph  of  the  initial  zero  flow  yields  a  graph 
with  no  negative  cvcles.  The  main  observation  of  Onaga  is  that  if  the  augmentation  is  done  using 
the  most-efficient  flow-generating  cycle,  then  all  flow-generating  cycles  in  the  residual  graph  of  the 
resulting  generalized  flow  pass  through  the  source.  Onaga’s  algorithm  iteratively  augments  the 
generalized  flow  along  the  most  efficient  flow-generating  cycle  in  the  residual  graph.  This  algo¬ 
rithm  maintains  the  invariant  that  every  vertex  is  reachable  from  s  in  the  current  residual  graph, 
a  property  that  can  be  verified  by  induction  on  the  number  of  augmentations. 

Consider  the  special  case  of  a  network  with  two  distinguished  vertices  s,t  €  V  and  with  all 
gains  other  than  7 {t,s)  and  7 (s,<)  equal  to  one.  The  generalized  flow  problem  in  this  network  is 
equivalent  to  the  maximum  flow  problem,  and  Onaga’s  algorithm  specializes  to  the  Ford-Fulkerson 
maximum  flow  algorithm.  Consequently,  the  algorithm  does  not  run  in  finite  time. 

The  next  algorithm  we  describe  uses  a  maximum  flow  computation  as  a  subroutine.  To  describe 
this  algorithm,  we  need  to  introduce  a  relabeling  technique  of  Glover  and  Klingman  [38].  This 
technique  can  be  motivated  as  follows.  Recall  the  financial  analysis  interpretation  of  the  generalized 
flow  problem,  in  which  vertices  correspond  to  different  securities  or  currencies,  and  arcs  correspond 
to  possible  transactions.  Suppose  one  country  decides  to  change  the  unit  of  currency.  (For  example. 
Great  Britain  could  decide  to  introduce  the  penny  as  the  basic  currency  unit,  instead  of  the  pound, 
or  Italy  could  decide  to  erase  a  couple  of  0’s  at  the  end  of  its  bills.)  This  causes  an  appropriate 
update  of  the  exchange  rates.  Some  of  the  capacities  change  as  well  (a  million  C  limit  on  the  C 
-  DM  exchange  would  now  read  as  a  limit  of  100  million  pennies).  It  is  easy  to  see  that  such  a 
relabeling  defines  an  equivalent  problem.  Such  a  relabeling  can  be  used,  for  example,  to  normalize 
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To  formally  define  the  relabeled  problem ,  let  p(v)  £  R+  denote  the  number  of  old  units  corre¬ 
sponding  to  each  new  unit  at  vertex  v  £  V.  Given  a  function  p,  we  shall  refer  to  p(v)  as  the  label 
of  v. 

Definition:  For  a  function  p  :  V  • — ►  R+  and  a  network  N  =  (V,E, 7,u),  the  relabeled  network  is 
wrhere  the  relabeled  capacities  and  the  relabeled  gains  are  defined  by 

u^r,  w)  =  u(v,w)/p{v) 

ln(v,w)  =  * y{v,w)p(v)/p(w ). 

For  a  generalized  pscudoflow  g  and  a  labeling  p,  the  generalized  pseudoflow  corresponding  to 
g  in  the  relabeled  network  is  g^{v,w)  =  g(v,w)/ p(v),  the  relabeled  residual  capacity  is  defined 
by  ug^(v,w)  =  (u(v,w)  —  g(v,w))/p(v),  and  the  relabeled  excess  by  e5iM(r)  —eg(v)/p(c).  The 
corresponding  pseudoflows  have  the  same  residual  graph. 

Now  we  present  a  canonical  relabeling.  The  residual  graph  of  a  generalized  pseudoflow  g  can  be 
canonically  relabeled  if  every  vertex  v  £  V  is  reachable  from  the  source  and  every  flow-generating 
cycle  in  the  residual  graph  contains  the  source.  For  a  vertex  v  £  V,  the  canonical  label  p(v)  is 
defined  to  be  the  gain  of  a  highest-gain  simple  path  from  s  to  v  in  the  residual  graph.  That  is. 
one  new  unit  corresponds  to  the  amount  of  flow  that  can  reach  the  vertex  t>  if  one  old  unit  of  flow 
is  pushed  along  a  most-efficient  simple  path  in  the  residual  graph  from  s  to  v,  ignoring  capacity 
restrictions  along  the  path.  Observe  that  in  a  restricted  network,  the  highest-gain  paths  from  s  to 
each  other  vertex  can  be  found  using  any  single-source  shortest  path  algorithm. 

Theorem  6.2.1  After  a  canonical  relabeling,  the  following  properties  hold:  Every  arc  £  Eg 

such  that  w  s  has  7M( v,w )  <  1;  there  exists  a  path  from  s  to  every  other  vertex  r  in  the  residual 
graph  with  7w(r,  it)  =  1  for  all  arcs  (i\w)  on  the  path;  the  most  efficient  flow-generating  cycles  each 
consist  of  an  arc  (r,  s)  £  Eg  and  an  (s,r )-path  for  some  r  £  V ,  such  that  7M(i\tr)  =  1  along  the  path 
and  7^(r, s)  =  max^^n,  w)  :  (v,u>)  £  Eg). 

Next  we  describe  a  simple  finite  algorithm,  due  to  Truemper  [101].  The  algorithm,  described  in 
Figure  6-1,  is  a  refinement  of  Onaga’s  algorithm  in  which  augmentation  along  all  of  the  maximum 
gain  cycles  is  done  simultaneously.  The  algorithm  maintains  a  generalized  pseudoflow  g  whose 
residual  graph  has  every  vertex  reachable  from  the  source  and  has  all  flow-generating  cycles  passing 
through  the  source.  The  algorithm  first  canonically  relabels  the  residual  graph.  Now  the  highest- 
gain  residual  cycles  have  maximum  relabeled  gain  on  arcs  entering  the  source  and  have  a  relabeled 
gain  of  one  on  all  other  arcs.  Consider  the  subgraph  induced  by  arcs  with  relabeled  gain  of  one 
and  arcs  of  maximum  relabeled  gain  entering  the  source.  A  circulation  that  maximizes  the  sum  of 
the  flow  on  the  latter  arcs  can  be  found  by  a  maximum  flow  computation.  This  circulation  gives 
an  augmentation  of  the  current  generalized  flow  g.  After  such  an  augmentation,  all  gain  cycles  in 
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Step  1  Find  /i,  a  canonical  labeling  from  the  source. 

If  yM(v,tu)  <  1  on  every  arc  of  the  residual  graph  of  the  current  generalized  flow  g, 
then  STOP  (the  current  flow  is  optimal). 

Step  2  Let  a  =  ma *(„,«>)€£,  7M(v,u>).  Consider  the  network  G'  consisting  of  all  arcs  (u,u>) 
with  JulVftu)  =  1,  and  all  arcs  in  A  =  {(t>,  s)lyu(v,  s)  =  a}  and  their  opposites. 
Find  a  (standard)  circulation  /'  in  G'  that  maximizes  ,)€A  f'{v,  s),  the  flow 
into  s. 

Step  3  Let  g'  be  the  generalized  flow  corresponding  to  /',  i.e,  g'(v,w)  =  f'(v,w)  if  v  ^  s 
and  g'(s ,  v)  =  -y^v,  s)f'(v,  s). 

Update  the  current  solution  by  setting  g(v,  u>)  =  g( v,  w)  +  g'(v,w)fi(v)  V(v,  w)  £  E. 


Figure  6.1:  The  Inner  loop  of  Truemper’s  algorithm. 


the  residual  graph  pass  through  the  source,  and  the  maximum  gain  of  a  flow-generating  cycle  in 
the  residual  graph  is  decreased.  Therefore,  this  algorithm  runs  in  finite  time. 

Truemper’s  algorithm  is  in  some  sense  an  analog  of  Jewell's  [61]  minimum-cost  flow  algorithm, 
which  augments  the  flow  along  all  of  the  cheaper!  augmenting  paths  at  once  using  a  maximum 
flow  subroutine.  Using  the  same  network  as  Zadeh  [105]  used  for  the  minimum-cost  flow  problem. 
Ruhe  [90]  gave  an  example  on  which  Truemper’s  algorithm  takes  exponential  time. 

Goldberg,  Plotkin,  and  Tardos  [42]  gave  two  algorithms,  one  that  uses  a  minimum-cost  flow 
computation  as  a  subroutine,  and  one  that  builds  on  ideas  from  several  maximum  and  minimum- 
cost  flow  algorithms.  In  the  next  section  we  describe  the  first  algorithm  in  detail  and  very  briefly 
outline  the  second  one. 


6.3  Polynomial-Time  Combinatorial  Algorithms 


The  main  idea  of  the  first  polynomial-time  algorithm  of  Goldberg-Plotkin-Tardos  is  best  described 
by  contrasting  the  algorithm  with  Truemper’s.  In  each  iteration,  both  algorithms  solve  a  sim¬ 
pler  flow  problem,  and  interpret  the  result  as  an  augmentation  in  the  generalized  flow  network. 
Truemper’s  algorithm  is  slow  because  at  each  iteration  it  considers  only  arcs  with  unit  relabeled 
gain  and  some  of  the  axes  adjacent  to  the  source,  disregarding  the  rest  of  the  graph  completely. 
The  Goldberg-Plotkin-Tardos  algorithm,  which  we  shall  call  algorithm  MCF ,  considers  all  arcs.  It 
assigns  a  cost  c(u,tn)  =  —  log  y^(v,  iv)  to  each  arc  (v,w)  and  so'vec  the  resulting  minimum-cost 
circulation  problem  (disregarding  gains). 

The  interpretation  of  a  pseudoflow  /  is  a  generalized  pseudoflow  g,  such  that  g(v,  w)  =  f(v,w)  if 
f(v,w)  >  0  and  g(v,w)  =  -y^w,  v)f(w,  v)  otherwise.  In  Truemper’s  algorithm,  the  interpretation 
of  a  feasible  circulation  on  the  subnetwork  G'  is  a  feasible  generalized  flow  on  the  original  network. 
In  Algorithm  MCF,  however,  the  interpretation  of  a  minimum-cost  circulation  is  a  generalized 
pseudoflov:.  Arc:  of  the  flow  that  have  a  relabeled  gain  of  less  than  1  produce  vertices  with  deficits 
in  the  interpretation.  A  connection  between  a  pseudoflow  /  and  its  interpretation  is  given  by  the 
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Step  1  Find  g,  a  canonical  labeling  of  the  residual  graph  of  the  generalized  pseudoflow  g. 
If  Tv(v.tc)  <  1  on  every  arc  of  the  rec.dual  graph  and  Vv  €  (V  —  {s})  :  e3iM(v)  =  0, 
then  HALT  (the  current  flow  is  optimal). 

Step  2  Introduce  costs  c(v,w)  =  —  log7>J(t u>)  on  the  arcs  of  the  network.  Solve  the 
capacitated  transshipment  problem  in  the  residual  relabeled  network  of  G  with 
demands  d(v)  =  eg  M(v)  for  t>  €  V  —  {s}. 

Step  3  Let  g'  be  the  interpretation  of  /'.  Update  the  current  solution  by  setting  g(v,w)  = 
g{ v,  u>)  4-  g'(v,  w)fi(v)  V(u,  u)  €  E. 


Figure  6.2:  Inner  loop  of  Algorithm  MCF. 


following  lemma. 


Lemma  6.3.1  The  residual  graphs  of  a  pseudoflow  /  and  its  interpretation  g  as  a  generalized  pseud¬ 
oflow  are  the  same. 

The  first  iteration  of  the  algorithm  solves  a  minimum-cost  circulation  problem,  and  it  creates 
excess  at  the  source  and  deficits  at  various  other  vertices.  Each  subsequent  iteration  uses  some  of 
the  excess  at  the  source  to  balance  the  deficits  created  by  the  previous  iterations  by  solving  a  ca¬ 
pacitated  transshipment  problem.  More  precisely.  Algorithm  MCF,  shown  in  Figure  6.2,  maintains 
a  generalized  pseudoflow  g  in  the  original  (non- relabeled)  network,  such  that  the  excess  at  every 
vertex  other  than  the  source  is  non-positive.  The  algorithm  proceeds  in  iterations.  At  eacn  itera¬ 
tion  it  canonically  relabels  the  residual  graph,  solves  the  corresponding  capacitated  transshipment 
problem  in  the  relabeled  network,  and  interprets  the  result  as  a  generalized  augmentation. 

The  most  important  property  of  a  minimum-cost  pseudoflow’,  for  this  application,  is  that  its 
residual  graph  has  no  negative  cycles.  This  implies  (by  Lemma  6.3.1)  that  the  residual  graphs  of 
the  generalized  pseudoflows  produced  by  the  algorithm  in  Step  3  have  no  flow-generating  cycles. 
The  following  lemma  (analogous  to  Lemma  1.6.1)  implies  that  this  is  enough  to  ensure  that  the 
new  residual  graph  can  be  canonically  relabeled. 

Lemma  6.3.2  Let  g  be  a  generalized  pseudoflow  in  a  restricted  network.  If  the  residual  graph  of  Gg 
contains  no  flow-generating  cycles,  and  the  excess  at  every  vertex  other  than  s  is  non-positive,  then 
every  vertex  is  reachable  from  s  in  G3. 

The  relabeled  gain  factors  are  at  most  1  (except  for  arcs  entering  the  source  in  the  first  iteration ); 
therefore  the  flow  computation  creates  deficits,  but  no  excesses  at  vertices  other  than  the  source. 
Using  the  decomposition  of  pseudoflows  (Theorem  1.7.3),  one  can  prove  that  in  each  relabeled 
network  there  exists  a  flow  satisfying  all  the  deficits.  These  properties  are  summarized  in  the 
following  lemma. 

Lemma  6.3.3  For  a  generalized  pseudoflow  g  that  is  constructed  by  Step  3,  the  following  properties 
hold:  The  residual  graph  of  g  has  no  flow-generating  cycle;  canonical  relabeling  applies  to  the  residual 
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graph  of  y;  all  excesses,  except  the  excess  of  the  source,  are  non-positive;  the  pseudoflow  required  in 
Step  2  of  the  next  iteration  exists. 

The  algorithm  terminates  once  a  generalized  flow  has  been  found.  Throughout  the  computation 
the  residual  graph  contains  no  flow-generating  cycles;  therefore  (by  Theorem  1.6.2)  if  the  algorithm 
ever  finds  a  generalized  flow,  then  this  how  is  optimal. 

Next  we  bound  the  number  of  iterations  of  the  algorithm.  Consider  the  generalized  pseudoflow  g 
existing  at  the  beginning  of  an  iteration.  Because  there  are  deficits  but  no  flow-generating  cycles  in 
the  residual  graph,  the  current  excess  at  the  source  is  an  overestimate  of  the  value  of  the  maximum 
possible  excess.  It  is  easy  to  see  that,  after  a  canonical  relabeling,  the  sum  of  the  deficits  at  all  the 
vertices  other  than  the  source  is  a  lower  bound  on  the  amount  of  the  overestimate.  We  use  this 
value,  De/[y,/z)  =  23„^a(— ea,e(r))'  as  a  measure  of  the  proximity  of  a  generalized  pseudoflow  to  an 

optimal  generalized  flow.  It  is  not  too  hard  to  show  that  if  Def{g,p)  is  smaller  than  L-1  then  the 
algorithm  terminates  in  one  more  iteration.  (Recall  that  L  is  the  product  of  denominators  of  arc 
gains.) 

Theorem  6.3.4  [42]  If  De/(y,/0  <  L-1  before  Step  2  of  an  iteration,  then  the  algorithm  produces  a 
generalized  flow  in  Step  3. 

An  important  observation  is  that  the  labels  p  are  monotonically  decreasing  during  the  algo¬ 
rithm.  The  decrease  in  the  labels  is  related  to  the  price  function  associated  with  the  minimum-cost 
pseudoflow'  computation.  Let  p'  denote  the  optimal  price  function  associated  with  the  pseudoflow 
f  found  in  Step  2.  Assume,  without  loss  of  generality,  that  p'{s)  =  0. 

Lemma  6.3.5  For  each  vertex  v,  the  canonical  relabeling  in  Step  1  of  the  next  iteration  decreases  the 
label  p(v)  by  at  least  a  factor  of  2~P'^VK 

The  key  idea  of  the  analysis  is  to  distinguish  two  cases:  Case  1,  in  which  the  pseudoflow  f  is 
along  “cheap”  paths,  (i.e.,  p'{v)  is  small,  say  p'lv )  <  log  1.5,  for  every  v  £  V);  and  Case  2,  in  which 
there  exists  a  vertex  v  such  that  p'(v)  is  “large”  (>  log  1.5).  In  the  first  case,  Def(f.p)  decreases 
significantly;  in  the  second  case,  Lemma  6.3.5  implies  that  at  least  one  of  the  labels  decreases 
significantly.  The  label  of  a  vertex  v  is  the  gain  of  the  current  most  efficient  path  from  s  to  r.  This 
limits  the  number  of  times  Case  2  can  occur.  Theorem  6.3.4  can  be  used  to  guarantee  that  Case  1 
cannot  occur  too  often.  The  following  lemma  provides  a  tool  for  estimating  the  total  deficit  created 
when  interpreting  a  minimum-cost  pseudoflow  as  a  generalized  pseudoflow. 

Lemma  6.3.6  Let  /'  be  a  pseudoflow  along  a  simple  path  P  from  ,s  to  some  other  vertex  v  that 
satisfies  one  unit  of  deficit  at  v.  Let  g'  be  the  interpretation  of  /'  as  a  generalized  pseudoflow.  Assume 
that  all  relabeled  gains  along  the  path  p  are  at  most  1,  and  denote  them  by  "n,...  7;  Then  after 
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augmenting  by  g' ,  the  unit  of  deficit  at  v  is  replaced  by  deficits  on  vertices  along  P  that  sum  up  to  at 

most  (rii<«<Jfc  7»)-1  ~  1- 

Proof:  The  deficit  created  at  the  ith  vertex  of  the  path  is  (1  -  7,),  for  i  =  1  ,...,Jt.  Using  the 
assumption  that  the  gain  factors  along  the  path  are  at  most  1,  the  sum  of  the  deficits  can  be 
bounded  by 


Eo -•*><£ 


1  -  7. 


•=1 


n;=i7,  nf=i7. 


- 1. 


This  lemma  can  be  used  to  bound  the  value  of  Def(g.p)  after  an  application  of  Step  3.  Let 
p'  be  an  optimal  price  function  associated  with  the  pseudoflow  /',  such  that  p'(s)  =  0.  Let  3  = 
maxt,ev  p'(u)  be  the  maximum  price.  Using  the  pseudoflow  decomposition  theorem  (Theorem  1.7.2). 
one  can  show  that  a  minimum-cost  pseudoflow  f  can  be  decomposed  into  flows  along  paths  from 
the  source  s  to  the  other  vertices  and  cycles  in  the  residual  graph  of  g ,  such  that  the  cost  o(  the 
cycles  is  zero  and  the  cost  of  each  path  ending  at  a  vertex  t»  is  at  most  p'(  v).  To  bound  the  amount 
of  deficit  created  by  the  interpretation,  we  shall  consider  the  interpretation  of  the  flow  on  the  cycles 
and  the  paths  one-by-one.  Since  the  cycles  consists  of  zero-cost  arcs  only,  the  interpretation  of  the 
flow  along  these  cycles  does  not  create  deficits.  When  interpreting  flow  along  a  path  from  ,«  to  v. 
the  deficit  at  v  is  replaced  by  deficits  that  sum  to  at  most  2P'^  -  1  <  23  -  1  times  the  deficit  at  v 
satisfied  by  this  path.  This  proves  the  following  lemma. 

Lemma  6.3.7  The  value  of  Def(g,p )  after  an  application  of  Step  3  can  be  bounded  by  23  -  1  times 
its  value  before  the  step.  In  particular,  if  <  log  1.5  for  every  vertex  v,  then  Def(g.p)  decreases 
by  a  factor  of  2. 

The  remain. ng  difficulty  is  the  fact  that  the  function  Def(g,fi )  can  increase,  both  when  Case  2 
applies  in  Step  3  and  due  to  the  relabeling  in  Step  1  (changing  the  currency  unit  from  C  to  penny 
increases  a  5  C  deficit  to  500  pennies).  The  increase  in  Def[g,p)  can  be  related  to  the  decrease 
in  the  labels,  however.  The  deficit  at  a  vertex  increases  in  Step  1  by  a  factor  of  a  if  and  only  if 
the  label  of  this  vertex  decreases  by  the  same  factor.  The  increase  in  Def(g.g)  during  Step  3  is 
bounded  by  23 ,  where  0  =  ma xp'(v),  by  Lemma  6.3.7.  By  Lemma  6.3.5,  this  means  that  /i(r)  for 
some  vertex  v  decreases  by  at  least  0  during  Step  1  of  the  next  iteration.  Hence,  in  both  cases,  if 
Def(g,p)  increases  by  a  during  a  step,  then  there  exists  a  vertex  v  for  which  p(v)  decreases  by  at 
least  a  during  the  next  execution  of  Step  1.  Let  T  <  Bn  denote  an  upper  bound  on  the  gain  of  any 
simple  path.  The  label  of  a  vertex  v  is  the  gain  of  a  simple  path  from  s  to  r  in  the  residual  graph. 
Therefore  the  labels  are  at  least  71-1  and  at  most  T.  and  the  label  of  a  vertex  cannot  decrease  by 
more  than  a  factor  of  T 2  during  the  algorithm.  This  gives  the  following  lemma. 
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Lemma  6.3.8  The  growth  of  the  function  Def(g,fi )  throughout  an  execution  of  the  algorithm  is 
bounded  by  a  factor  of  T =  B°(n2\ 


Theorem  6.3.9  [42]  The  algorithm  terminates  in  0(n2  log  B)  iterations. 

Proof:  The  label  of  any  vertex  v  is  the  gain  of  a  simple  path  from  s  to  v,  and  it  is  monotonically 
decreasing  during  the  algorithm.  Therefore,  it  cannot  decrease  by  a  constant  factor  more  than 
O(logT)  times  for  any  given  vertex,  and  Case  2  cannot  occur  more  than  O(nlogT)  times.  When 
Case  1  applies,  the  value  of  Def(g,/i )  decreases  by  a  factor  of  2.  The  value  of  Def(g,fi)  is  at  most 
O(nBT)  after  the  first  iteration;  and,  by  Theorem  6.3.4,  the  algorithm  terminates  when  this  value 
decreases  below  L~l .  Lemma  6.3.8  limits  the  increase  of  Def(g,g)  to  T°^  during  the  algorithm. 
Hence,  Case  1  cannot  occur  more  than  0(log(nBT)  +  nlogT  +  log  L)  =  0[n2  log  B )  times.  | 

To  get  a  bound  on  the  running  time,  one  has  to  decide  which  algorithm  to  use  for  computing 
the  minimum-cost  pseudoflow  in  Step  2.  The  best  choice  is  Orlin's  0(rr.(m  +  «logn)logn)  time 
algorithm,  discussed  in  Chapter  5. 

Theorem  6.3.10  [42]  Algorithm  MCF  solves  the  generalized  flow  problem  using  at  most  0(n2m(m  -f 
n  log ti )  log n  log B )  arithmetic  operations  on  numbers  whose  size  is  bounded  by  O(m\og  B). 

We  conclude  this  section  with  a  brief  discussion  of  the  other  algorithm  of  [42].  The  algorithm  is 
based  on  ideas  from  two  flow  algorithms:  the  cycle-canceling  algorithm  of  Goldberg  and  Tarjan  [45] 
described  in  Chapter  4,  and  the  fat-path  maximum  flow  algorithm  of  Edmonds  and  Karp  [21],  which 
finds  the  maximum  flow  in  a  network  by  repeatedly  augmenting  along  a  highest-capacity  residual 
path. 

Each  iteration  of  the  algorithm  starts  with  cycle-canceling.  Canceling  a  flow-generating  cycle  in 
a  generalized  flow  network  creates  an  excess  at  some  vertex  of  the  cycle.  If  we  cancel  cycles  other 
than  the  most  efficient  ones,  the  residual  graph  of  the  resulting  generalized  pseudoflow  will  have 
flow-generating  cycles  that  do  not  contain  the  source.  The  algorithm  cancels  flow-generating  cycles 
by  an  adaptation  of  the  Goldberg-Tarjan  algorithm.  Using  the  dynamic  tree  data  structure,  this 
phase  can  be  implemented  to  run  in  0(n2m log n  log  B)  time.  The  resulting  excesses  at  vertices 
other  than  the  source  are  moved  to  the  source  along  augmenting  paths  in  the  second  phase  of  the 
algorithm. 

Consider  a  generalized  pseudoflow  that  has  non-negative  excesses  and  whose  residual  graph  has 
no  flow-generating  cycles.  The  key  idea  of  the  second  phase  is  to  search  for  augmenting  paths  from 
vertices  with  excess  to  the  source  that  result  in  a  significant  increase  in  the  excess  of  the  source. 
The  algorithm  maintains  a  scaling  parameter  A  such  that  the  maximum  excess  at  the  source  is  at 
most  2mA  more  than  the  current  excess.  It  looks  for  an  augmenting  path  that  can  increase  the 
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excess  of  the  source  by  at  least  A.  If  the  residual  graph  of  the  current  generalized  pseudoflow  has 
no  flow-generating  cycles,  then  one  can  find  a  sequence  of  such  paths  such  that,  after  augmenting 
the  generalized  flow  along  these  paths,  the  maximum  excess  at  the  source  is  at  most  mA  more 
than  the  current  excess.  This  phase  can  be  implemented  in  0(m(m  +  nlogn))  time  by  a  sequence 
of  O(m)  single-source  shortest  path  computations  on  graphs  with  arcs  of  non-negative  cost.  Now 
A  is  divided  by  two  and  a  new  iteration  is  started.  After  O(m\og  B)  iterations,  when  the  excess 
at  the  source  is  very  close  to  the  optimum  value,  Truemper’s  algorithm  can  be  applied  to  bring  all 
of  the  remaining  excesses  to  the  source. 

Maggs  (personal  communication,  1989)  has  observed  that  this  algorithm  can  be  improved 
through  better  balancing  of  the  two  phases.  Dividing  A  by  a  factor  of  2a  after  each  iteration 
decreases  the  number  of  iterations  by  a  factor  of  a  and  increases  the  time  required  for  the  second 
phase  by  a  factor  of  2° .  The  parameter  a  is  chosen  so  that  the  time  required  for  the  two  phases  is 

balanced.  The  resulting  algorithm  runs  in  O  (n2m2  log2  i?log((ns ^gn^HTogiogfi)  lime- 
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